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ABSTRACT 


This  report  describes  a  methodology  and  means  for  accurately 
determining  confidence  limits  for  the  reliability  of  large  digital 
logic  networks,  without  exhaustively  exercising  all  possible  input 
sequences  or  simulating  all  logic  faults  in  the  network.  The  re¬ 
port  is  divided  into  two  parts.  In  the  first  part,  some  mathemat¬ 
ical  models  to  estimate  the  reliability  of  digital  circuits  are 
presented.  A  heuristic  method  of  assigning  weights  to  faults  de¬ 
pending  on  their  “importance" 'in  a  given  circuit  is  described. 

The  models  presented  can  be  used  to:  Ca)  predict  the  reliability 
of  a  circuit,  (b)  evaluate  test  sequences  and  (c)  develop  more 
accurate  reliability  models  of  (redundant)  fault  tolerant  compu¬ 


ters.  The  second  part  of  the  report  deals  with  the  statistical 


methods  of  estimating  confidence  limits. 
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Chapter  I 
INTRODUCTION 


With  the  increased  complexity  of  current  digital  systems,  reliability 
considerations  have  become  increasingly  important.  Being  physical  devices, 
digital  circuits  are  subject  to  failure.  Although  current  technologies  em¬ 
ployed  to  construct  digital  systems  are  more  reliable  than  earlier  tech¬ 
nologies,  to  a  great  extent  the  resulting  decrease  in  the  failure  rate  of 
individual  components  has  been  offset  by  the  increased  complexity  of  today's 
circuits.  This  is  one  of  the  main  reasons  for  the  increased  interest,  both 
in  industry  and  academia,  in  the  subjects  of  maintenance  and  reliability  of 
digital  systems. 

A  fault  or  failure  of  a  digital  network  can  be  defined  as  a  physical 
defect  of  one  or  more  components  in  the  network  which  causes  it  to  behave 
differently  from  the  original  systems  (Su  [74]).  In  digital  systems,  typ¬ 
ical  maintenance  goals  deal  with  the  rapid  detection,  location  and  repair 
of  any  system  faults.  In  many  digital  systems  involving  real-time  process¬ 
es,  such  as  telephone  switching  networks  and  aircraft  or  spacecraft  flight 
controls,  it  is  desirable  to  continuously  monitor,  exercise  and  test  the 
system  in  order  to  determine  whether  the  system  is  performing  as  desired. 
Such  monitoring  may  enable  automatic  detection  of  failures  via  periodic 
testing  or  through  the  use  of  codes  and  checking  circuits  (e.g.  self-test¬ 
ing  and  self-checking  circuits)  or  may  enable  continuous  operation  under 
failure  (i.e.,  fault  tolerence)  and  automatic  repair  via  switching  networks 
(e.g.  stand-by  spares). 


One  way  to  determine  whether  a  fault  exists  in  a  circuit  is  to  exercise 
the  network  against  all  possible  input  sequences  and  compare  the  results  with 
expected  output.  A  second  method  is  to  insert  faults,  physical  or  modelled, 


into  a  copy  of  the  network,  then  generate  tests  that  detect  the  presence  of 
these  faults.  For  a  complex  digital  system  such  as  an  airborne  computer, 
it  is  not  only  impossible  to  apply  all  possible  input  sequences  but  also 
impractical  to  model  all  possible  faults.  A  more  practical  and  accepted 
approach  is  to  generate  input  patterns  that  detect  a  certain  percent  of 
faults.  The  generation  and  evaluation  of  such  test  sequences  have  been 
the  subjects  of  much  investigation  (see  for  example  Friedman  [71]  or  Breuer 
[76]). 

As  the  average  number  of  IC's  on  a  board  increases,  the  difficulty  of 
adequately  testing  the  board  increases,  perhaps  exponentially.  It  can  be 
shown  that  with  normal  rejection  rates  of  incoming  IC's,  boards  of  50  to  100 
IC's  will  always  contain  at  least  one  bad  IC  prior  to  testing  (Fike  [72]). 

As  things  stand  now,  checking  out  a  $10  microprocessor  chip  may  re¬ 
quire  an  investment  by  the  user  of  upwards  of  1000  times  that  amount  in 
test  gear  (e.g.  Fairchild's  Sentry  VII)  -  and  he  still  may  not  know  if  the 
chip  will  do  the  job  it  is  supposed  to  do  (Vodovoz  [75]).  Both  user  and 
manufacturer  come  face  to  face  with  the  problems  of  seeing  if  the  chips 
work  and  at  this  moment,  no  one  is  fully  satisfied  with  present  chip  check¬ 
out  techniques.  The  hardest  hit,  though,  is  the  end  user  who  measures  his 
chip  needs  in  hundreds,  and  who  can  not  afford  the  sophisticated  test  equip¬ 
ment  that  can  give  him  better  answers  than  those  he  gets  with  his  own  home 
brewed  test  methods. 

Two  truisms  must  be  recognized  at  this  point: 

1.  Testing  does  not  add  to  product  quality,  it  merely  evaluates  the  quality 
already  present. 

2.  The  test  function  normally  requires  a  larger  expenditure  for  equipment 
than  needed  for  any  other  part  of  the  production  organization. 
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Another  way  of  expressing  the  first  statement  is  to  say  that  the  qual¬ 
ity  must  be  built  in  through  use  of  good  components  and  good  workmanship. 
The  second  statement  indicates  that  the  question  of  the  amount  of  testing 
to  be  performed  must  be  very  carefully  considered  and  answered.  The  amount 
of  testing  which  should  be  performed  is  that  minimum  which  clearly  demon¬ 
strates  that  the  product  performs  or  fails  to  perform  as  specified.  This 
rather  vague  definition  must  be  expanded  greatly  so  that  all  responsible 
personnel  have  a  reasonably  accurate  concept  as  to  what  constitutes  a  'pro¬ 
per  minimum  test'  and  also  what  constitutes  the  economics  of  testing.  One 
such  attempt  appeared  in  Watkins  [70]. 

Most  users  demand  that  the  test  equipment  detect  troubles  which  are 
not  present  at  the  time  of  testing  but  which  may  occur  at  some  future  time. 
To  achieve  such  tight  tolerences,  one  has  to  test  his  system  on  a  regular 
basis.  The  cost  of  such  testing  can  be  prohibitively  expensive.  So,  we 
need  a  procedure  which  enables  the  end  user  to  evaluate  the  trade-offs 
between  the  frequency  of  testing  (consequently,  the  cost  of  testing)  and 
the  "confidence"  he  can  have  in  his  system.  This  is  the  subject  of  our 
Investigations. 

This  report  describes  a  methodology  and  means  for  accurately  deter¬ 
mining  confidence  limits  for  the  reliability  of  large  digital  logic  net¬ 
works,  without  exhaustively  exercising  all  possible  input  sequences  or 
simulating  all  logic  faults  in  the  network.  This  report  is  divided  into 
two  parts.  In  the  first  part,  some  mathematical  models  to  estimate  the 
reliability  of  digital  circuits  are  presented.  A  heuristic  method  of  as¬ 
signing  weights  to  faults  depending  on  their  "importance"  in  a  given  cir¬ 
cuit  is  described.  The  models  presented  can  be  used  to 
1.  predict  the  reliability  of  a  circuit 
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2.  evaluate  test  sequences 

3.  develop  more  accurate  reliability  models  of  (redundant)  fault  tolerent 
computers. 

The  second  part  deals  with  the  statistical  methods  of  estimating  confi¬ 
dence  limits. 
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Chapter  II 


TEST  SEQUENCES 

Test  generation  is  the  process  of  finding  the  set  of  input  patterns 
for  a  digital  circuit  which  will  either  verify  that  the  circuit  is  operat¬ 
ing  correctly  or  else  provide  some  information  on  the  nature  of  the  failure. 

The  set  of  such  input  patterns  is  often  designated  as  a  test  sequence  (or 
simply,  a  test).  It  is  desirable  that  the  test 

1.  be  reasonable  in  size 

2.  be  produced  at  reasonable  expense,  and  yet, 

3.  detect  a  maximum  number  of  faults. 

Many  algorithms  have  been  proposed  for  generating  tests  for  digital 
circuits.  Some  of  the  procedures  make  use  of  an  algebraic  description  of 
the  circuit  under  consideration.  Others  directly  utilize  the  gate  level 
circuit  topology  and  functional  description.  In  this  section  we  will  sur¬ 
vey  briefly  some  of  these  algorithms.  In  almost  all  these  studies,  a  com¬ 
mon  assumption  is  that  a  logic  element  can  only  fail  by  sticking  at  zero 
or  by  sticking  at  one. 

The  test  generation  procedures  that  have  been  evolved  so  far  can  be 
broadly  classified  as  either  deterministic  or  probabilistic.  Examples  of 
deterministic  techniques  include  the  Boolean  difference  method  (Sellars  [68]),  the 
one  dimensional  path  sensitization  method  (Armstrong  [66])  and  D-algorithm 
(Roth  [66]).  All  of  these  methods  offer  high  fault  coverage  (i.e.,  each 

test  pattern  detects  a  loarge  number  of  faults)  and  reasonable  size  test 
sequences.  But,  they  are  computationally  complex  and  limited  to  special 
classes  of  circuits.  For  example,  the  D-algorithm  uses  cubical  algebra  and 
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requires  extensive  programming  to  automate  the  procedure.  Moreover,  the  D- 
algorithm  is  originally  designed  for  combinational  circuits.  Extensions 
to  these  methods  to  handle  sequential  circuits  have  evolved  in  recent 
years,  but  the  complexity  prohibits  their  use. 

In  the  random  technique,  a  candidate  test  is  chosen  using  a  random 
number  generator.  If  the  candidate  test  detects  new  faults  not  detected 
by  previous  tests,  it  is  added  to  the  test  set;  otherwise  it  is  discarded. 

The  figure  of  merit  of  a  candidate  test  can  be  defined  as  the  ratio  of  the 
number  of  new  faults  detected  to  the  number  of  faults  to  be  detected.  A 
candidate  test  may  be  discarded  if  its  figure  of  merit  is  less  than  a  spec¬ 
ified  threshold  value. 

There  are  two  major  categories  of  Automatic  Test  Equipments  (ATE),  call¬ 
ed  stored  program  ATE  and  comparison  ATE.  Stored  program  testers  usually 
contain  a  mini  computer  and  back-up  storage  as  disk  and  test  sequences  stor¬ 
ed  vector  by  vector  or  as  a  high  level  program  interpreted  by  the  computer.  The 
stored  program  ATE  typically  also  stores  the  responses  and  a  fault  diction¬ 
ary,  which  are  usually  generated  by  simulation.  The  actual  test  sequences 
can  be  obtained  using  either  a  deterministic  or  random  procedure. 

Comparison  ATE  employs  pseudo -random  patterns  as  test  vectors.  Here 
two  circuits,  the  Unit  Under  Test  (UUT)  and  a  known  good  copy  of  the  circuit 
(denoted  by  C*)  are  inserted  into  the  ATE  (Figure  1) . 
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Figure  1.  Comparison  ATE 

The  pattern  generator  applies  several  million  pseudo  random  patterns  to  both 
£ 

the  UUT  and  C  .  The  outputs  of  the  circuits  are  then  compared  by  the  discrep¬ 
ancy  detector.  A  mismatch  indicates  a  fault  in  UUT.  The  advantages  of  com¬ 
parison  type  ATE  are  1.  very  fast  generation  of  test  patterns,  2.  absence 
of  an  expensive  computer  and  storage.  However,  a  known  good  copy  of  the  cir¬ 
cuit  is  needed. 

Although  random  test  patterns  can  be  generated  very  inexpensively,  this 
technique  becomes  progressively  inefficient  when  attempting  to  detect  more 
deeply  embedded  faults  (Parker  [75a]).  Adaptive  (Parker  [75a])  and  Weighted 
(Schnurmann  [75])  random  test  generations  have  been  developed  to  improve  the 
efficiency  of  random  test  generation. 

For  a  complex  digital  network,  one  is  usually  satisfied  with  tests  that 
detect  a  certain  percent  of  the  modelled  faults  (say  95%) .  As  this  number 
approaches  100%,  the  number  of  input  patterns  required  to  detect  these  faults 
approaches  a  very  large  number.  So,  most  end  users  are  satisfied  with  tests 
that  check  a  percentage  of  faults. 
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Let  G  denote  the  set  of  all  modelled  faults  and  let  Q  be  the  set  of 
faults  an  user  want  to  detect  (Q<G) .  Typically,  if  an  end  user  is  satis¬ 
fied  with  tests  that  detect  95%  of  the  modelled  faults,  then  95%  of  G  consti¬ 
tutes  Q.  A  test  sequence  T  is  then  generated  to  detect  the  faults  in  Q 
and  the  circuit  is  exercised  against  all  the  input  patterns  in  T  periodical¬ 
ly  to  assure  reliable  operation  of  the  circuit . 

Let  Tl,  T2,...,Tn  be  a  number  of  test  sequences  detecting  faults  in  Ql, 
Q2,...,Qn  where  each  Qi  is  95%  (or  any  given  percent)  of  G.  The  actual  set 
of  faults  detected  by  each  of  these  tests  Ti  can  be  determined  before  hand 
by  probing.  By  using  a  fault  dictionary,  a  fault  may  be  identified  with¬ 
in  its  equivalent  class.  The  set  of  physical  faults  associated  with  this 
equivalent  class  may  be  distributed  over  many  components,  such  as  IC's.  In 
order  to  determine  the  exact  location  of  a  physical  fault,  probing  techniques 
are  required.  It  should  be  clear  that  by  interchanging  these  test  Ti  peri¬ 
odically,  say  Tl  during  the  first  testing,  T2  when  the  circuit  is  tested 
second  time  and  so  on,  the  reliability  associated  with  the  circuit  can  be 
enhanced  appreciably  (Figure  2) . 


Figure  2.  Sets  of  Faults 
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PROBABILITY  OF  DETECTING  A  FAULT 

Parker  [75b]  has  defined  signal  probability  as  follows: 

Definition:  The  probability  of  a  signal  denoted  as 

a  -  P(A-l)  (1) 

for  a  signal  A  is  the  probability  that  signal  A  equals  1.  Similarly, 
the  probability  that  the  signal  equals  0  is  given  by 

P(A-O)  -  1  -  P(A-l)  =  1  -  a  (2) 

Boolean  operations  NOT,  AND,  OR  can  be  applied  on  probabilities. 

(a)  Boolean  Negation  (NOT)  corresponds  to  the  probability  expression 

b  -  1  -  a  (3) 

(b)  Boolean  AND  of  two  independent  signals  A  and  B  in  the  expression  C  * 
A.B  corresponds  to  the  probability  expression 

c  *  a  •  b  (4) 

(c)  Boolean  OR  of  two  independent  signals  A  and  B  in  the  expression  C  * 

A  +  B  corresponds  to  the  probability  expression 

c  «  a  +  b  (5) 

Let  f (xl,x2, . . . ,xn)  denote  the  function  realized  by  a  combinational 
circuit.  A  logical  fault  a  changes  the  function  realized  to  f  (xl,x2,... 
xn) .  Using  the  theory  of  boolean  difference  methoc  (Sellers  [68]),  the 
fault  a  is  detected  by  an  input  for  which  f  +  f  is  1.  Using  the  probabil¬ 
ity  expression  for  NOT,  AND,  OR  operations,  the  probability  P(f  +  f  -1)  can 


be  evaluated  as  a  function  of  the  input  signal  probabilities.  This  is  the 
probability  that  fault  a  will  be  detected.  It  is  often  convenient  to  assign 


equal  probabilities  to  input  signals.  This  is  known  as  bundling  (Parker 
[75c]). 

To  see  how  the  probability  of  detecting  a  fault  can  be  used  in  random 
testing  let  us  consider  the  circuit  in  Figure  3. 


f  -  ABC  +  ABC  +  ABC  +  ABC  (6) 

Figure  3  An  Example  -  Probability  of  detecting  a  fault. 

Let  a  be  the  fault  that  A  is  stuck  at  1.  Then 


f0  (A,B,C) 

»  B  +  C 

(7) 

f  +  f 

a 

■  BC  +  ABC 

(8) 

P(f  +  f 

■1)  =  b  +  c  -  ac  -  2bc  +  abc 

(9) 

Assuming  equal  probabilities  to  a,  b,  and  c,  i.e., 
a  ■  b  *  c  *  x 

the  probability  of  detecting  fault  a  is  given  by 

P(f  +  f  *  1)  ■  2x  -  3x^  +  x^  (10) 

a 

Since  x  is  a  probability,  the  above  function  is  plotted  for  values  of  x 
between  0  and  1. 
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p 


0.4 


x 


This  function  P(f  +  f  ■  1)  in  terms  of  the  input  signal  probabilities 
can  be  used  to  evaluate  the  values  cf  the  input  signal  probabilities  for 
any  desired  probability  of  detecting  the  fault  a.  Theses  values  can  be 
used  to  control  the  random  input  patterns  generated  leading  to  adoptive 
test  generation.  For  example,  from  the  above  plot,  the  input  signal  prob¬ 
abilities  (i.e.,  x)  must  be  set  to  0.4  in  order  to  achieve  a  probability 
of  detecting  the  fault  of  0.4. 

DEPENDENCY  OF  RANDOM  TEST  GENERATION  ON  THE  LOGIC  DEPTH 

The  Monte  Carlo  method  was  used  extensively  for  testing  the  circuits  of 
ILLIAC  IV  (Moreno  [72],  Agrawal  [72]).  Some  of  the  interesting  results 
obtained  during  the  testing  are  given  here.  It  was  found  that  a  circuit 
with  437  lines  was  completely  tested  (all  s.a.O  and  s.a.l  faults  with  sin¬ 
gle  fault  assumption)  with  56  random  input  patterns.  Another  circuit  with 
89  lines  required  210  patterns  for  complete  testing  while  for  a  third  cir¬ 
cuit  with  115  lines,  the  test  generation  was  not  complete  even  after  2000 
patterns.  These  results  lead  to  the  conclusion  that  the  computation  time 
does  not  depend  only  on  the  number  of  lines  or  faults.  Agrawal  [75a]  showed 
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that  the  number  of  random  inputs  required  for  a  complete  test  of  a  circuit 
depends  on  the  logic  depth,  i.e.,  the  number  of  levels  in  the  circuit.  In 
Agrawal  [75b] ,  an  expression  for  the  probability  of  sensitizing  a  path  upto 
the  primary  output  through  L  levels,  P(L)  is  derived.  Then,  the  probability 
of  sensitizing  a  path  through  L  levels  by  at  least  one  out  of  M  independent 
patterns  P(L,M)  is  given  by 

P(L,M)  -  1  -  [1  -  P(L)JM  (11) 

Solving  the  above  equation  for  M, 


M  _  l2£lL.7-P.a,M)  ] 

log[l  -  P(L) ] 


(12) 


This  equation  can  be  used  to  estimate  the  number  of  random  input  patterns 
required  to  detect  all  faults  in  a  circuit. 

Experimental  results  have  shown  that  equation  (12)  gives  a  good  estimate 
of  the  number  of  random  input  patterns  required  for  a  complete  check  of  com¬ 
binational  circuits  (Agrawal  [75b]).  Although  these  results  are  derived  for 
simple  tree  structures  consisting  of  n  input  NAND  gates,  the  model  can  be 
applied  to  any  combinational  circuit  since  an  equivalent  NAND  trees  can  be 
constructed  for  the  given  circuit  very  easily. 


CONFIDENCE  LEVEL 

In  Agrawal  [75b],  the  term  confidence  is  used  to  mean  the  probability 
of  sensitizing  a  path  through  L  levels  by  at  least  one  of  the  M  random  in¬ 
put  patterns  in  a  test  set.  Another  way  of  describing  this  meaning  is  to 
say  that  the  confidence  of  a  test  sequence  is  the  probability  of  detecting 
a  fault  since  a  fault  is  detected  if  a  path  from  the  site  of  the  fault  to 
a  primary  output  can  be  sensitized. 
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In  another  attempt  to  define  confidence  level,  Shedletsky  [77]  derived 
an  expression  for  latency  intervals  in  a  circuit.  Error  latency  can  be  de¬ 
fined  as  the  delay  between  the  occurrence  of  a  fault  and  the  first  error  in 
the  output.  Latency  interval  of  a  circuit  is  the  maximum  of  the  minimum 
number  of  input  patterns  necessary  to  achieve  a  given  probability  of  detect¬ 
ing  a  fault.  Shedletsky  notes  that  the  required  length  of  a  random  test 
(i.e.,  the  number  of  input  sequences)  to  achieve  a  given  confidence  level 
is  equal  to  the  latency  Interval  of  the  circuit.  A  more  detailed  discussion 
of  error  latency  is  included  in  a  later  section.  Thus,  in  this  definition, 
confidence  level  is  related  to  the  number  of  test  patterns  in  that  test  set. 

In  both  of  the  above  attempts,  confidence  level  is  defined  in  terms  of 
a  test  sequences.  However,  a  more  useful  definition  of  confidence  level 
should  relate  to  the  reliability  of  the  circuit  itself.  A  meaningful  def¬ 
inition  should  include  mean  time  between  failures  (mtbf)  of  the  circuit  due 
to  logical  and  physical  faults.  Such  a  definition  can  be  used  to  schedule 
periodic  check-outs  in  order  to  achieve  a  desired  confidence  in  the  circuit. 
In  the  next  section  two  definitions  of  reliability  are  given.  The  authors 
believe  that  these  definitions  are  more  useful. 


RELIABILITY  ASSOCIATED  WITH  A  TEST  SEQUENCE 

Here  two  definitions  of  reliability  are  presented.  The  first  defini¬ 
tion  accounts  for  the  faults  that  are  detected  by  the  test  sequence  under 
consideration.  This  definition  is  then  extended  to  Include  modelled  faults 
that  are  not  detected  by  the  test  sequence  and  even  faults  that  can  not  be 
modelled . 

Definition  1:  The  reliability  of  a  circuit  associated  with  a  test  se¬ 
quence  T  which  detects  the  set  of  faults  Q  -  (gl,g2 , . . . ,gn)  is  a  function 
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of  time  and  is  given  by 


n  [ (1-R.  /c  fi'dt)] 
i-1  1  to 


(13) 


assuming  that  the  circuit  has  passed  the  test  T  at  time  tO  and  where 

fi  is  the  probability  density  function  of  the  fault  gi  with  a  mean  of 

1/X± 

and  ki  is  a  constant  which  describes  the  dependence  of  reliability 
on  the  fault  gi,  0  <  ki  <  1. 

The  term  (1  -  fi-dt)  in  the  above  equation  gives  the  probability 

that  the  fault  gi  does  not  exist  at  time  t  since  /t  fi*dt  is  the  cumulative 

to 

probability  that  fault  gi  is  present  at  time  t  assuming  that  gi  did  not 
exist  at  time  tO. 

Typically,  all  ki’s  are  set  equal  to  1,  however,  an  user  can  bias  the 
reliability  in  favor  of  some  of  the  faults  depending  on  their  importance  in 
his  circuit.  These  constants  are  in  line  with  the  weighted  random  test  pat¬ 
tern  generation  (Schnurmann  [75])  where  a  weight  to  each  signal  is  assigned 
according  to  its  importance  in  the  circuit. 

Equation  (13)  describes  some  kind  of  decay  function  for  R.J,  and  Figure  5 
shows  the  general  shape  of  such  a  curve.  By  applying  the  test  T  periodically, 
the  reliability  can  be  restored  to  the  maximum  (provided  the  circuit  pass¬ 
es  the  test  T  each  time).  In  such  cases,  the  curve  has  a  sawtooth  form 
(Figure  6) .  The  frequency  at  which  the  circuit  must  be  tested  depends  on 
the  minimum  reliability  Rj,. 
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Figure  5  Reliability  Associated  with  a  Test. 


Time 

Figure  6  Periodic  Testing  Restores  Reliability  to 
its  Maximum 
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The  above  definition  takes  into  account  only  the  faults  that  are  de¬ 
tected  by  the  test  T.  But,  intuitively,  it  should  be  obvious  that  the  re¬ 
liability  of  the  circuit  should  be  lower  than  that  predicted  by  the  above 
equation  due  to  the  faults  that  are  not  detected  by  test  T.  So,  a  second 
definition  of  the  term  reliability  is  needed. 

Let  ■  ^®n+l’®n+2’ . gn)  t*ie  set  m°delled  faults  that  are 

not  detected  by  test  sequence  T.  Let  NM  =  (gnmi»gnm2 . ,gnTTl1 )  be  the 

set  of  unmodelled  faults.  Certain  physical  defects  like  shorts  can  be 
modelled  as  logical  faults  while  some  physical  defects  like  changes  in  vol¬ 
tage  and  loading  can  not  be  modelled  as  logical  faults.  Defects  that  can 
not  be  modelled  may  effect  the  performance  of  the  digital  system. 

Definition  2:  The  reliability  of  a  circuit  associated  with  a  test  sequence 
T  is  given  by 

n  nm£  £ 

r»  [  n  (i-ki  /c  fi-dt)][  n  (l  -k,  /c  f  dt)][  n  u- /  f.-dt)] 

i*l  to  j»n+l  tj  J  £»nml  to 

where  fi,  f j ,  and  fl  are  the  probability  density  functions, 
ki  and  kj  are  weighting  constants  as  before  and 
tO  is  the  time  the  test  T  was  applied  last, 
tj  is  the  time  when  a  test  checked  for  fault  gj  and 
t  is  the  time  when  the  circuit  was  designated  operative. 

The  first  term  in  equation  (14)  is  due  to  the  faults  (gl,g2, . . . ,gn) 
that  are  detected  by  test  T.  The  second  term  is  due  to  the  modelled  faults 
that  are  not  detected  by  test  T  and  the  third  term  is  due  to  the  unmodelled 
faults.  The  third  term  may  be  dropped  if  the  unmodelled  faults  are  not 
critical  or  if  the  failure  rates  of  such  faults  are  not  readily  available. 

The  general  form  of  the  reliability  curve  given  by  equation  (14)  is 
shown  in  figure  7.  When  the  reliability  reaches  a  value  which  is  not  ac- 


(14) 


16 


cepCable,  the  circuit  must  be  replaced. 


time 

The  first  definition  may  be  used  when  an  user  is  interested  only  in 
the  faults  that  are  detected  by  his  test  set.  But,  if  more  accurate  re¬ 
liability  measures  are  desired,  the  second  definition  should  be  used. 


AN  EXAMPLE 


The  reliability  calculations  are  illustrated  using  the  following 
circuit. 


Figure  8.  An  Example 


This  circuit  has  24  stuck  faults  under  the  single  fault  assumption.  This  number 
is  reduced  to  10  faults  on  input  signals  and  4  faults  on  branches  at  fanout 
points  (g  and  h  stuck  faults)  using  fault  collapsing.  The  table  given  below 
contains  6  input  patterns  that  detect  all  faults  in  the  circuit. 


Inputs _ _ Faults  Detected 


ABODE 

1 

110  10 

a  b  c  e.  f  h.  i.  j.  k.  r.  m. 

0-  o,  0,  1,  O,  1,  1,1,  1,  1,  1 

2 

0  10  10 

a,  c.  f,  g-  J  m 

1,  1,  1,1,  O,  O 

3 

10  0  10 

b.  c-  f-  g-  J  m 

1,  1,  1,1,  o,  0 

4 

0  *  1  1  0 
x  0 

c  d  h  i  k  m 

O  y  O  y  O  y  O  )  Of  O 

5 

x  0 

10  0 

0  x 

d,  el  g  i,  J.  k,  m, 

1,  ,Bo,  1,  1,  1,  1 

6 

x  x  1  0  1 

e  k  m 
o,  o,  o 

Table  -  1 

Let  fa  *  i  e  ^  be  the  probability  density  function  (pdf)  of  the 
o 

fault  a  stuck  at  0.  For  the  sake  of  this  example,  we  will  assume  that  all 
faults  have  identical  pdf's.  We  will  assign  a  1  to  the  weighting  constants. 
Test  set  T  *  (1,3,5)  detects  75%  of  the  faults  in  the  circuit. 


Qt  «  (a0,al,b0,bl,c0,cl,dl,el,f0,fl,g0,gl,hl,il, jO, jl,kl,m0,ml) 
Using  the  first  definition, 

o  n  / 1  ,t.  -Ar.  .  -19At 

R.J,  -  (1  -  /  Ae  dr)  -  e 


If  the  mean  time  between  failures  (1/X)  has  a  value  of  5  X  10  hours,  then 

„  -19x2xlO~5t 

®T  6 

For  a  minimum  reliability  of  0.9,  the  time  between  testings  can  be  calculated: 
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t  ■  lnCR.jJ/C-SS  x  10  )  ■  277  hours 

Thus,  the  test  T  should  be  applied  every  277  hours  in  order  to  maintain  a 
reliability  of  0.9. 

While  using  the  second  definition,  we  will  ignore  faults  that  can  not 
be  modelled.  The  set  of  faults  that  are  not  detected  by  test  set  T  is 

0^  -  (d0,e0,h0,i0,k0) . 

Assuming  that  all  faults  were  tested  at  time  t*0, 


_-Xr 


R-,  *  [  II  (1  -  /  Xe  dr)  ]  [  H  (1  -  Xe  Aidr)] 
T  i=l  o  1=20 


When  X  *  2  x  10  ,  for  a  reliability  of  0.9,  the  time  between  tests  is  cal¬ 

culated  to  be  t  a  219.5  hours. 

Table  (2)  lists  few  other  values  for  time  between  test  for  various 
reliabilities. 


TABLE  2 


50.000  hrs. 


2  x  10 
2  x  10“ 
2  x  10“ 
2  x  10“ 


1.14  x  10 
1.14  x  10 
1.14  x  10 


ti  with  1st  def. 


757  hrs 


ti  with  2nd  def. 


559  hrs 


Chapter  III 

ASSIGNMENT  OF  WEIGHTS  TO  FAULTS 


In  a  microprocessor  system  or  in  any  digital  system,  certain  faults 
are  more  crucial  to  the  operation  of  the  system  than  others.  For  example, 
in  most  IC's,  the  enable  and  power  lines  are  the  most  critical  faults 
bn  these  are  fatal  to  the  system.  It  may  be  useful  to  assign 
weights  to  the  signals  according  to  their  importance  in  a  circuit  so 
that  better  reliability  measures  can  be  obtained.  It  is  very  difficult  to 
come  up  with  a  procedure  which  assigns  absolute  values  to  weights  since  the 
importance  of  a  signal  not  only  depends  on  its  function  in  a  circuit,  but 
the  relative  importance  of  signals  may  be  biased  by  user  views.  Here  we 
attempt  to  outline  a  heuristic  procedure  that  enables  an  user  to  assign 
weights  to  various  faults. 

Faults  in  a  digital  system  may  be  classified  into  four  groups. 

1.  Very  Important  Faults:  These  faults  are  very  critical  to  the  operation 
of  the  system.  It  may  be  useful  to  include  faults  whose  relative  importance 
is  not  entirely  clear.  This  leads  to  a  conservative  estimate  of  the  circuit 
reliability,  but  the  weights  may  be  changed  when  the  significance  of  faults 
becomes  discernible.  We  can  assign  a  weight  of  1  to  faults  in  this  class. 

2.  Important  Faults:  There  may  be  some  faults  whose  presence  impairs  the 
operation  of  the  circuit  but  does  not  create  a  hazardous  outcome.  These 
faults  may  become  critical  to  the  circuit  operation  or  prolonged  existence. 

For  example,  a  stuck  input  in  a  shift  register  does  not  effect  the  output 
immediately,  but  if  the  iault  remains  it  may  lead  to  an  erroneous  results. 

Such  faults  are  assigned  a  weight  of  0.75. 

3.  Unimportant  Faults:  Faults  in  this  class  may  not  be  critical  to  the 
operation  of  the  circuit  at  all,  even  in  continued  existence.  For  example, 

in  most  computers,  a  fault  in  a  clock  circuit  which  changes  the  cycle  time  slightly 
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is  not  critical  to  the  operation  of  the  system.  A  weight  of  0.5  is  assign¬ 
ed  to  faults  in  this  category. 

4.  Don't  Care  Fault:  Some  faults  may  not  effect  the  operation  of  the  cir¬ 
cuit,  especially  if  redundant  logic  is  used.  These  faults  can  be  ignored  in 
some  dircumstances,  and  so  a  weight  of  0.0  is  assigned  to  such  faults. 

At  this  point  we  would  like  to  draw  an  analogy  from  an  automobile. 
Examples  of  very  important  faults  include  faults  in  the  circuit  that  con¬ 
trols  the  brakes,  faults  in  the  ignition  circuit.  A  fault  in  spedometer 
or  a  gasoline  indicator  can  be  critical  in  prolonged  existence.  A  faulty 
spare  tire  can  also  be  classified  as  an  important  fault.  A  failure  in  the 
heater  or  air  conditioner  operation  can  be  classified  as  unimportant  faults. 

A  defect  in  accessories,  such  as  the  radio  or  clock,  can  be  considered  as 
a  don't  care  faults. 

BOUNDS  ON  TIME  BETWEEN  TESTS 

In  this  section  the  upper  and  lower  bounds  on  the  time  between  tests 
(TBT)  are  derived.  The  reliability  of  a  circuit  associated  with  a  test  set 
S  drops  below  the  permissible  minimum  if  the  circuit  is  not  exercised  against 
the  test  set  S  at  least  once  every  TBTm^y  for  that  test.  Application  of  a 
test  more  frequently  than  that  indicated  by  TBTm;tn  may  not  be  necessary  to 
maintain  the  desired  reliability. 

Let  Sl,S2,...,Sn  be  a  number  of  test  sequences  such  that 

S-US-US-U. . .US  *  G 

1  2  3  n  (15) 

SL  *  Sj  if  iH 

Let  us  consider  a  test  sequence  S.  Let  g  be  the  most  important  fault  among 

the  faults  detected  by  S.  Let  t  denote  the  time  between  test  for  test  set  S. 

s 


Upper  Bound  on  tg:  The  maximum  value  of  tg  satisfies  the  following  equation 

1  -  kS  CS  fS'dt  3  ^  (16) 

where  k  is  the  weight  assigned  to  fault  g,  fg  is  the  pdf  associated  with 
© 

g  and  is  the  minimum  acceptable  reliability. 

Lower  Bound  on  gs :  Let  Qg  *  (gl,g2, . . .gm)  be  the  set  of  faults  detected  by 
the  test  set  S.  Assuming  that  all  the  faults  are  equally  important  in  the 
circuit,  we  can  assign  a  weight  of  1  to  all  faults  in  Qs.  Then,  the  minimum 
value  of  ts  satisfies 

ill  C1  -  £Sl'dt)  ■  <17> 

The  actual  value  of  the  time  between  tests  for  the  test  set  S,  tg, 
depends  on  the  real  values  of  the  weights  k^'s.  If  reasonable  weights  can 
not  be  assigned  to  faults  in  a  digital  system,  TBTmin  should  be  used  for 
scheduling  the  check-outs  of  the  circuits. 

FAILURE  RATES  OF  VARIOUS  FAULTS 

Sometimes  it  may  not  be  practical  to  calculate  the  failure  rates  of 
each  and  every  fault  in  a  digital  system.  In  such  cases,  it  is  desirable 
to  classify  the  faults  into  various  groups  depending  on  either  the  proxim¬ 
ity  of  the  nature  of  the  failure  mechanisms  or  on  the  closeness  of  the 
failure  rates.  Once  similar  faults  are  grouped  together,  identical  fail¬ 
ure  rates  for  faults  in  a  group  may  be  assumed.  Since  the  grouping  of 
faults  into  various  categories  impacts  the  accuracy  of  the  reliability  mea¬ 
sures,  the  classification  should  be  conducted  very  carefully. 

Let  be  the  failure  rates  respectively  of  the  faults  gl, 

g2,...,gn,  the  modelled  faults  that  are  detected  by  the  test  T.  Let 
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Xn+^,Xn+2, . . . Xr  be  the  respective  failure  rates  of  the  modelled  faults 

gn+i*gn+2»  •  •  •  >Sn  tliat  are  not  detected  by  the  test  T  while  >*nTn? » •  • « 

X  .  are  the  failure  rates  of  unmodelled  faults, 
max 

Assinning  exponential  failure  rates,  the  sum  of  these  failure  rates 

gives  the  failure  rate  of  the  digital  system  under  consideration, 
n  r  nm£ 

Xckt  *  l  +  l  Xj  +  l  X£  (18) 

i*l  j*n+l  £*nml 


Since  it  is  not  possible  to  find  the  individual  failure  rates  of  unmodel¬ 
led  faults  (or  their  nature) ,  it  is  reasonable  to  conglomerate  the  failure 

rates  of  unmodelled  faults  into  a  single  term  Xnm. 
n  i 

Xckt  "  I  U  +  l  Xj  +  Xnm  (19) 

i-1  j-n+1 


The  fault  rates  of  modelled  faults  may  similary  be  bundled  into  one  or  more 
groups  depending  on  the  similarity  of  the  faults  and  the  required  accuracy 
of  the  modal.  We  will  illustrate  bundling  of  faults  for  LSI  circuits. 

For  LSI  circuits  it  is  observed  (Tees  [71],  Kasouf  [78])  that  the 
failure  rates  of  the  pins  (primary  inputs  and  outputs)  are  higher  than  the 
failure  rates  of  the  gates.  This  leads  to  a  classification  of  faults  into 
3  groups  -  faults  on  pins,  faults  on  gates  and  unmodelled  faults.  The  fail¬ 
ure  rate  can  be  written  as 


Lct  *  A  .  +x_  +x 

LSI  pins  gates  nm 


(20) 


where  X  is  the  total  failure  rate  of  faults  on  pins, 

Xgateg  is  the  failure  rate  due  to  faults  on  gates. 

Equation  (20)  is  similar  to  the  equations  in  Tees  [71]  and  Kasouf  [78] 
which  are  of  the  general  form 

lLSI  *  VP  +  V0  +  C3  (21) 
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where  P  is  Che  number  of  pins  and  G  is  the  number  of  gates  in  a  LSI  circuit. 

Constant  Cl,  C2,  and  C3  are  evaluated  from  the  data  supplied  by  the  vendor 

from  his  own  failure  analysis  and  user  feedback.  In  our  case,  the  model 

is  similar,  but  the  constants  reflect  the  failure  rates  of  modelled  logic 

faults  on  pins  and  gates.  If  the  values  of  X  and  \  in  equation 

pins  gauCS 

(20)  can  be  obtained,  failure  rates  of  individual  faults  can  be  calculated 
assuming  identical  fault  rates  for  faults  in  a  group.  These  failure  rates 
can  then  be  used  to  calculate  reliability  of  the  circuit  and  the  frequency 
of  testing  to  achieve  a  desired  reliability. 

Available  data  about  failures  in  semiconductor  devices  is  typically 
of  the  form,  n  failures  in  m  hours  of  operation.  This  type  of  data  enables 
us  to  estimate  the  failure  rate  of  the  circuit  as  a  whole.  However, 

if,  in  addition  to  observing  that  a  failure  has  occurred,  the  nature  of 
the  fault  -  a  fault  on  a  pin,  a  fault  on  a  gate  or  else  -  is  discovered, 
we  would  have  data  to  estimate  *pin  and  ^gates  in  equation  (20) .  If  the 
bundling  results  in  a  different  model,  the  failure  rates  due  to  faults  in 
each  group  can  similary  be  estimated. 


Chapter  IV 


EXTENSIONS  AND  SOME  OTHER  APPLICATIONS 


Although  equation  (20)  is  derived  using  the  LSI  failure  mechanism, 
this  model  (equation  (19))  can  be  used  for  any  circuit.  For  a  complex  circuit 
using  a  large  number  of  circuits,  two  extensions  of  the  model  are  possible. 

1. 

XCkt  -  l  XICi  +  |  Xwire  j  +  Xm  «« 

where 

X^  ^  is  the  failure  rate  of  IC  i  considered  as  one  piece. 

*wire  j  is  the  failure  rate  of  faults  on  the  interconnecting  wire 
j  which  depends  on  the  wire  length  etc., 

X  _  is  the  failure  rate  of  all  the  other  faults. 

rim 


x  ,  ■  x  +  x  .  +  x  +x 

Ckt  pins  wires  gates  nm 


(23) 


where 


Xpins  is  the  combined  failure  rate  of  faults  on  the  pins  of 
all  the  IC's. 

Xwires  is  the  combined  failure  rate  due  to  the  faults  on 
interconnecting  wires. 

Xgates  is  the  combined  failure  rate  of  the  faults  on  all 
gates . 

Xnm  is  the  failure  rate  due  to  all  the  other  faults. 


As  noted  in  the  introduction,  the  model  can  be  used  to  evaluate  test 
sequences.  Recently,  a  statistical  method  for  test  sequence  evaluation  is 
developed  (case  [16]).  Generation  and  evaluation  of  test  sequences  have  been 
the  subjects  of  much  Investigation  in  the  past.  But  here,  test  sequence 
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evaluation  is  a  by-product. 


In  the  previous  sections ,  we  have  shown  how  to  calculate  the  reliability 
associated  with  a  test  sequence.  This  reliability  can  be  used  as  a  measure 
to  rank  test  sequences:  a  better  test  will  have  higher  reliability  and/or 
lower  frequency  of  testing  for  a  given  confidence  level. 

One  other  possible  application  is  to  use  our  model  in  deriving  reliability 
equations  for  redundant  fault-tolerant  systems .  There  have  been  many  mathemati¬ 
cal  models  developed  for  redundant  systems  like  triple  modular  redundant  (TMR) , 
hybrid  redundant,  and  stand-by  spare.  (See  for  example  Bouricius  [71],  Mathur 
[71]).  In  almost  all  these  studies,  a  simple  exponential  function  is  used  to  rep¬ 
resent  reliability  of  constituent  systems.  If  the  model  presented  in  this  report  is 
used  instead  of  simple  exponential  distribution,  more  accurate  reliability  pre¬ 
dictions  can  be  derived. 
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CONFIDENCE  INTERVAL  FOR  THE  RELIABILITY  R. 

A 


As  the  reliability  measure  R^  defined  earlier  is  decreasing  as  time  in¬ 
creases,  we  would  like  to  maintain  the  reliability  at  some  level  by  applying 
the  test  sequence  A  periodically.  In  order  to  find  the  associated  confidence 
limit  to  maintain  RA  at  some  given  level  and  finally  to  present  the  time  inter¬ 
val  which  is  required  to  check  the  given  circuit  periodically,  let's  define 
some  terminologies  first. 

Definition 

Confidence  level  100  6%  (0  <_  S  <_  1)  of  the  reliability  is  the  percentage 
value  of  the  probability  such  that  we  are  at  least  100  BZ  sure  that  the  reli¬ 
ability  is  contained  in  some  interval  (R^.R^  which  is  called  the  100  BZ  con¬ 
fidence  interval  of  R^  while  and  Ru  are  called  lower  and  upper  confidence 
limit  of  RA  respectively,  i.e.  Prob.  [R^  <  R  <  Ru]  »S. 

Remark:  Frequently  one  is  interested  in  one-sided  confidence  interval  and  the 
definition  given  above  can  be  modified  accordingly  i.e.  if  one  Is  interested 
in  lower  confidence  limit  only  (as  is  the  case  with  us)  then  Prob.  [P^<  R]  ■  6 
is  the  associated  probability  statement  and  100  6%  is  called  the  confidence 
level  of  Ra  and  R^  is  called  the  (lower)  confidence  limit  of  RA  etc. 

If  Identical  failure  rates  are  assumed  among  each  of  QA>  QA  and  NM,  the 
reliability  RA  will  eventually  be  a  function  of  time  t  with  parameter 
X  »  X^  +  X0+X^  where  X^:  failure  rate  for  modeled,  detected  by  test  seq.  A. 

XjJ  failure  rate  for  modeled,  undetected  by  test  seq.  A. 

X^:  failure  rate  for  modeled,  unmodeled  faults. 

i.e.  Ra  •  R^(t;  X^^,!^)  ■  Rj(t;X)  where  X  *  X^  +  +  X^.  (23) 

Further  assume  exponential  failure  model  for  each  of  QA  ’  ^A  an<*  c^e  Para~ 


meters  X ^  and  X^  respectively  and  assume  failures  in  each  set  Q^»Qj  and 

occur  independently  of  each  other. 

Here,  one  performs  an  experiment  and  observes  a  Poisson  process  in  which 


NM 


r  failures  are  observed  in  time  t  ■  n  t  ,  where  t  is  the  time  of  termination 

r  o  o 


of  the  life  test  and  n  is  the  number  of  items  in  the  sample  subjected  to  life 
test  (Type  I  Censoring  See  Mann  [75])  and  r  *  ri+r2+r3’  where  r^,^  and  r^  are 
the  number  of  faults  in  Q.,Qt  and  NM  respectively  observed  in  time  interval  t  . 

A  A  O 

Let  K.  be  the  random  variable  of  which  r.  is  the  realization,  i  ■  1,2,3  and 
let  K  ■  +  K^.  Since  K^'s  follow  Poisson  distribution  with  parameter 

Xi  tr  (see  Mann  [75]),  K  also  follows  Poisson  with  parameter  Xt  by  the  inde¬ 
pendency  assumption  of  K^'s  for  i  •  1,  2  and  3. 


Let  0  *  n  t  /r  ■  t  ^r  and  using  the  identity  Prob  (K  r)»Prob[x  (2r+2)>2X  t^] , 


where  x  (v)  is  the  random  variable  representing  chi-square  distribution  with 

egrees  of  f 
X£_0  (2r+2) 


v  degrees  of  freedom,  we  can  show  100  (l-a)%  upper  confidence  limit  Xy  for  X 


is 


,  where  x  (v)  means  q-th  quantile  of  chisquare  distribution 
2r  9  q 


with  v  degrees  of  freedom.  Here,  one  can  obtain  the  lower  confidence  limit 


to  maintain  at  confidence  level  100(l-a)%,  since  is  a  decreasing  function 


of  y  .  i.e. 


®L  *  *  ®2  t"’ 


*i-a(2r+2) 


(24) 


2r  0 


where  R2  is  given  in  (23) 

* 

Further,  to  find  the  time  period  t  required  to  periodic  checking  procedure 
to  maintain  the  reliability  at  some  given  level  Rffl^n,  we  can  simply  use  (24) 
i.e.,  t*  is  the  value  of  t  which  satisfies  the  equation 


min 


t; 


xI-a(2r+2> 

A 

2r  0 


(25) 
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* 

i.e.  if  we  check  the  given  circuit  at  t  time  period,  we  can  at  least  be 

O'" 

100  (l-o)%  sure  that  the  reliability  R.  is  at  least  R  . 

A  min 


Note:  Since  we  do  not  know  the  true  value  of  R^,  which  is  a  function  of 
unknown  X^'s,  it  is  more  realistic  to  estimate  the  lower  bound  of 
R<j,  instead  of  "mean"  value  of  R^,  and  use  the  time  interval  to  main¬ 
tain  the  lower  bound  of  R^.  In  other  words,  results  in  this  sec¬ 
tion  is  more  conservative  (more  frequent  checking)  compared  to  the 
case  in  the  examples  in  the  appendix  at  the  end  of  Chapter  II. 
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Chapter  V 


Point  and  Interval  Estimation  for  Life  Testing 
Procedure  of  LSI/MSI  Reliability  Models  -  System  Point  of  View 

In  theoretical  studies  of  equipment  reliability,  one  is  often  concern¬ 
ed  with  systems  consisting  of  many  components,  each  subject  to  an  individual 
pattern  of  malfunction  and  replacement,  and  all  parts  together  making  up  the 
failure  pattern  of  the  equipment  as  a  whole. 

Drenick  [1960]  showed  that  a  complex  piece  of  equipment,  after  an  extend¬ 
ed  period  of  operation,  will  tend  to  exhibit  a  failure  pattern  with  an  expoen- 
tial  distribution  for  the  time  between  failure  (inter-arrival  times)  and  that 
the  time  up  to  the  first  failure  is  also  nearly  exponentially  distributed. 

Hence  LSI/MSI  reliability  model  can  be  described  by  expoential  density  function 
f(t)=*  (l/0)exp(-t/0)  (1) 

where  9  is  the  average  life  time  and  t  is  the  time  to  first  failure 

of  LSI/MSI. 

One  is  often  interested  in  estimating  0  or,  if  mission  time  tffl  is  speci¬ 
fied,  in  estimating  reliability  of  an  LSI /MSI  component  R(tm)  *  1  -  F(tm)  where 

F(t)  is  the  cumulative  distribution  function  for  f(t).  i.e.  R(t  )  is  the  pro- 

m 

bability  of  no  fault  occurs  in  the  given  LSI/MSI  until  time  tm> 

In  the  following  notes,  the  estimation  procedure  of  0  and  R(-)  will  be  explain¬ 
ed  and  demonstrated  through  examples  for  the  various  sampling  situations.  More¬ 
over  the  procedure  of  computing  the  time  period  to  maintain  the  given  LSI/MSI 
component  at  the  specified  reliability  level  (lower  confidence  limit  for  R(O) 
is  illustrated. 

1.  Testing  without  replacement 

In  the  life  testing  procedure  of  LSI/MSI,  a  faulty  LSI/MSI  is  not  replac¬ 
ed  until  the  testing  is  terminated. 


(a)  Estimation  for  Type  II  censoring 

Suppose  that  a  sample  of  size  n  LSI (MSI)  with  failure  time  distribution 
given  by  (1)  has  been  subjected  to  life  testing  and  that  the  test  is  termina¬ 
ted  at  the  time  of  the  rc^  failure,  r  £  n.  More  generally,  one  might  assume 
that  a  sample  of  size  n  has  been  randomly  selected  from  a  one-parameter  expo¬ 
nential  population,  and  that  experimentation  is  terminated  at  the  time  that 
rck  observation  becomes  available. 


The  joint  density  function  of  the  ordered  observations  X^,, 
<_  »  i*l,...,r-l  and  r*l,2,...,n  is  given  by 


X(r)  ’ 


Si)  ’  ‘  "  ’X(r)  (X1’  *  *  *  ,Xr}*  (n-rjle1 


iI1xi+(n-r)xr 


O  £  Xi  .  .  ,<^  Xr 


Clearly,  if  the  right-hand  side  of  this  is  maximized  with  respect  to  9, 
one  obtains,  as  the  maximum  -  likelihood  estimator  (MLE)  of  9, 


i£lx(i)-Kn-r)X(r) 


Epstein (1953-1954)  showed  that  9  is  unbiased,  sufficient  and  complete  for 
9  and  because  9  is  based  on  a  complete  sufficient  statistic,  it  is,  by  the 
Lehmann-Scheffe* -Blackwell  theorem,  unique  minimum-variance  estimator  (UMVUE) 
among  unbiased  estimator  of  0. 

Also  Var(9)  ■  9^/r. 

If  we  let  *  (n-i+l)(X^  -  i*l,...r,  has  an  Independent 

expoential  distribution  with  scale  parameter  equal  to  1  and  it  is  well  known 
that  2S^,  i»l,...r,  has  an  independent  chi-square  distribution  with  2  degrees 


of  freedom. 


r 

Also  one  can  easily  infer  that  E  2S.  ■  2r0/0  (being  the  sum  of  r  indepen- 

i=l  1 

dent  chi-squares)  has  a  chi-square  distribution  with  2r  degrees  of  freedom.  A 
lower  confidence  bound  on  0  at  confidence  level  1-a  is  therefore  given  by 

»  2  2 

2r9/x^  a(2r),  where  x^Ck)  is  the  100  th  percentile  of  a  chi-square  distribution 

*  9  *  2 

with  k  degrees  of  freedom.  The  interval  [2r9/x“_a/2  (2r)  ,  (2r)  ]  is  a 

two-sided  confidence  interval  for  9  at  confidence  level  1-a.  To  obtain  such 

intervals  one  simply  substitutes  calculated  values  for  0  and  tabulated  values 

of  chi-square  percentile  (see  example  below) .  The  minimum  variance  unbiased 

point  estimator  of  R(t  )  for  r=n  is  [See  Pugh  (1963)  and  Basu  (1964)] 

m 

R*(t  )  «f(l  -  ^r)n_1  ,  n0  >  t  , 
m  ]  n0  m 

L  0  *  n9  -  v 

A  lower  confidence  bound  for  R(t  )  at  level  1-a  can  be  obtained  by  substi- 

m 

*  2 

tuting  2r8/x^_a(2r)  for  S  in  the  expression  exp(-tm^g).  Two-sided  confidence 
intervals  are  obtained  similarly. 

These  results  are  summarized  in  Table  I  below. 

The  time  period  t*  to  maintain  the  reliability  at  a  given  level,  say  6, 


can  be  obtained  by  equating  0  with  the  lower  confidence  limit  for  R(tm) 

solving  for  t  . 

m 

Table  I 
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Example  1. 

Let  n*24  LSI(MSI)'s  are  put  in  life  testing  and  the  test  was  terminated 
at  r*3rd  failure. 

If  the  observed  exponential  failure  times  in  hours  from  the  above  censor¬ 
ed  samples  of  size  24  are  6200,9200  and  16900  then  9**[6.2  +  9.2  +  22(16.9)]* 
10^/3  =*  129,100  hours.  95%  two-sided  confidence  interval  (a*0.05)  is  [53605.5, 

626040.6].  If  mission  time  t  *  5,000  hours,  lower  confidence  limit  for  R(t  ) 

m  m 

at  level  .95  (a=.05)  with  t  *5,000  hours  is  exp  (-5000/ [^v^rTiT  1  )g .92  i.e. 
we  can  be  at  least  95%  sure  that  the  reliability  is  at  least  .92.  To  obtain 
the  time  period  to  maintain  the  reliability  at  least  8*. 95  level,  we  solve 

.-Wex^ioo, 

.95  *  exp  [  /  — 12 — 59 — ^  and  get 

'ft 

tffl  =  3155.8  i.e.  we  might  have  to  check  the  LSI  after  using  3155  hours  period¬ 
ically  to  be  95%  sure  that  the  reliability  is  at  least  .95. 


(b)  Estimation  for  Type  I  censoring 

If  a  life  test  is  terminated  at  a  specified  time  tQ,  the  joint  p.d.f. 
of  x^  , . . .  is  given  by 


(1) 


...,x(r)(xr...,xr)  -  r  exP  *i  +  (n-r)  tQ] } , 


0<X.<. . .  <X  <t  . 
—  1—  —  r —  o 


The  maximum  likelihood  estimator  of  9  is  easily  seen  to  be 
r 

i  ■  ,  tM. 

r  ’ 

This  estimator  is  biased  for  small  samples  but  has  all  the  desirable 
asymptotic  properties  associated  with  the  MLE  of  9  under  Type  II  censoring. 
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To  obtain  a  confidence  interval  for  9,  one  may  consider  the  probability 
that  the  number  of  failures,  K,  in  a  sample  of  size  n  is  equal  to  r  at  time  tD 

P(K~r)  =  (n-r) !r!  [R(t°> ^ [1'R(to) ]'• 

Then  a  conservative  (because  K  is  discrete  and  the  failure  times  are  not 
used)  lower  confidence  bound  6  at  level  1-a  for  R(t0)  is  given  by  the  solution 
of 


„  _  r(n+l) 


o  r(n-r)T(r+l) 


n-r+1,,  N r , 
x  (1-x)  dx 


From  the  relationship  between  the  binomial  and  the  beta  distributions, 
it  can  be  seen  that  6  is  equal  to 

1  -  V._  (r+1,  n-r)  =  V^Cn-r,  r+1) ,  the  lOOath  percentile  of  the  beta 

distribution  with  parameters  n-r  and  r+1.  From  this,  one  can  determine,  from 
R(t0)=exp(~to/9) ,  a  lower  confidence  bound  for  9  at  level  1-a  as  t0/ln(l/<5). 


Table  II 


parameter 


point  estimate 


interval  estimate 


i£l  X(i)  +  (n~r)to 


,  r  *  0 


to/ln(l/6) 

(lower  confidence  bound) 


R(t0) 


R(t;.0) 


6  =  V  (n-r,  r+1) 
a 

(lower) 


Example  2. 

Sample  size  n*15  LSI's  were  put  in  life  testing  and  the  test  was  termin¬ 
ated  at  to*240  days.  The  number  of  failures  up  to  time  t  were  r*2;  then  a 
90  %  lower  confidence  bound  for  R(tQ)  is  .6827  and  a  90  to  lower  confidence 


bound  for  9  is  240/ln(l/ .6827)*  630  days 


f  , . . . ,x,  . (x  , . . .x  )=(nA)rexp(-nAx  ) ,0<x  . .<x  , 

X(1)  V  ;  1  r  1  r 

where  A  =  1/0 . 

From  this  expression  one  can  obtain  the  m.L.E.  of  A  as  (r/X(r))/n 
and  of  0  as  (nX^rp/r.  For  this  model,  one  can  also  write  the  joint  p.d.f. 
of  x^  , » * .  as 

r  r 

(nA)  exp[-nA  E  (x.-x.  -)],  x  =0<x1£. . . <x 
i=-l  1 

In  other  words,  each  Si  *  X^-X^  i»l,...,r,  with  X=0,  has  an  inde¬ 
pendent  exponential  distribution  with  scale  parameter  (An)  Since  X^  * 

r 

E  Si,  the  density  of  X,  .  for  testing  with  replacement  is  given  by 
i-1  W 

/  ,^r  y  n_^ 

f  (xr)  =  -  f-  —  exp(-nAxr). 

X(r)  U  ij. 

Consequently,  2nAX(r)  has  a  chi-square  distribution  with  2r  degrees 
of  freedom,  and  confidence  bounds  for  A  ,  9  or  reliability  can  be  obtained 
as  functions  of  observed  values  X^r^  and  values  to  be  obtained  from  tables 
of  the  chi-square  distribution. 
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Table  III 


parameter 

point 

interval 

9 

o>  > 

If 

a 

/■“N 

ft 

fT 

2nX(r)  2nx(r) 

x2  a  X2a 

|_  l-f(2r)  f(2r) 

R(t) 

R(t;9) 

Mt; 

Xl-a(2r) 

(lower  limit) 

Example  3 

For  a  sample  (with  replacement)  of  size  5  LSI's  from  an  exponential 
population  for  which  the  first  observed  failure  time  is  5600  hours  and  the 
fourth  and  last  observed  time,  occuring  at  the  termination  of  the  life  test, 
is  20,000  hours,  the  MLE  for  9  is  5(20,000) /4**25, 000  hours  and  an  80  %  lower 
confidence  bound  for  9  is  10(20,000) /X2g0(8)=18, 100  hours. 

(b)  Estimation  for  Type  I  censoring. 

In  this  case,  one  observe  a  Poisson  process  in  which  r  failures  are 
observed  in  time  tr»nt0,  where  tQ  is  the  time  of  termination  of  the  life 
test  and  n  is  the  number  of  items  in  the  sample  subjected  to  life  test. 

For  this  model,  the  number  of  observed  failures,  K,  is  a  random  variate  with 

P(K«r)«(-iy)(\tr)rexp(-Atr)  with  A=l/9. 

The  M.L.E  X  of  A  is  r/tr  *  r/ntQ  and  the  MlE  9  of  9  is  A-^»nt0/r. 

then 

P(K<r)»  l  i  (Atr)kexp(-Atr) 
k-0 

-  zr 

*  f\tr  rT  exP('z)dz’ 
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’••***. 

An  incomplete  gamma  function  with  shape  parameter  r  +  1.  Therefore 
Prob(K<r)  *  Prob  [x2(2r+2)  >  2Xtr].  Thus  if  one  observes  r  failures,  then, 
with  probability  at  least  1-a,  2Ar9<X2  a(2r+2),  or  2nt0/x^_a (2r+2)  is  a  con¬ 
servative  lower  confidence  bound  for  9  at  confidence  level  1-a.  Similarly, 
it  can  be  shown  that 

Prob  (K>r)  =  Prob[X2(2r)  2Xtr]. 

Hence  a  conservative  upper  confidence  bound  for  9  at  level  1-a  is  given 
by  2r9/x^(2r),  and  a  conservative  two-sided  confidence  interval  at  level 
1-a  is 

2r9 _  2r9 

Xl-./2<2r+2>  ,  Xa/2(2r) 

“  m 

Cox  (1953)  discusses  the  interval 


[- 


2r9 


2r9 


Xl-a/2<2r+1> 


X;/2(2r+D 


which  he  states  is  slightly  narrower  than  that  given  above  but  sometimes 
has  a  true  confidence  coefficient  less  than  1-a. 


Table  IV 


parameter 

point 

interval 

9 

9-nt0/r 

2r9  ,  2r9 

_XL,/2<2r+2>  Xa/2(2t>_ 

R(t) 

R(t)  -  R(t;9) 

R(t.  -  ) 

X^_a(2r+2) 

(lower  limit) 

Note  that  11(^9)  is  an  increasing  function  of  9. 
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Example  4 

If  200  LSI's  are  life  tested  with  replacement  for  1,000  days  and  fail¬ 
ures  are  observed,  a  90%  lower  confidence  bound  for  the  mean  time  to  fail¬ 
ure  is  400  (1,000) /X2q(8)»400, 000/13. 36  *  29,900  days.  As  an  example  in 
real  life,  we  illustrate  the  estimation  procedure  using  the  data  from  G. 
Kasouf  and  S.  Mercurio  (1978). 

Certified  data  from  an  airborne  radar  processing  systems  (RPS)  operat¬ 
ing  in  a  simulated  airborne  inhabited  environment  is  used  to  calculate  the 
observed  LSI/MSI  failure  rate.  The  certified  data  was  accumulated  during 
the  RPS  1600-hour  test-analyze-and-f ix  program  followed  by  a  351-hour  fixed 
length  reliability  demonstration  test.  Analysis  of  data  shows  an  accumula¬ 


tion  of  3.9  and  4.6  million  operating  hours  for  LSI  and  MSI  circuits,  respec¬ 
tively.  One  LSI  and  no  MSI  failures  were  experienced.  In  this  case  tr* 

(3.9  +  4.6)  x  106  *  8.5  x  106  and  X«l/(8.5xl06)=0.118/106  hrs  and  r-1. 

60%  (a*. 4)  confidence  interval  for  X (failure  rate)  is 


X22(2)X  X28(4)r 

_  2  ’  2  \ 


[.026314  ,  .35341] 


95%  lower  confidence  bound  for  0  is 


X^95(4) 


10° 

2  x  0.118 
9.488 


1.78637  x  10°  hours 


To  obtain  the  time  period  to  maintain  the  LSI  at  B».98  level,  we  solve 


and  get 


.98  -  exp[-t/(l. 78637  x  106)] 


t*  ■  3 . 609  x  10^  hours 

m 


i.e.  we  may  have  to  check  the  LSI's  in  3.6  x  104  hours  period  to  be  95%  sure 


that  the  reliability  is  at  least  .98. 
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Appendix  A 


STATISTICAL  PROPERTIES  OF  QUASI-RANGE  IN  SMALL 
SAMPLES  FROM  A  GAMMA  DENSITY 
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Introduction 


Karl  Pearson  [  1920] ,  studied  combinations  of  order  statistics 
which  were  later  named  quasi-ranges  by  Mosteller  [1946] .  Quasi-ranges 
are  simple  to  obtain  and,  for  moderate  sample  sizes, yield  more  effi¬ 
cient  estimators  of  standard  deviation  than  do  sample  ranges.  In  cases 
where  abundant  data  are  available  and  where  the  cost  of  complicated 
data  reduction  far  outweighs  the  cost  of  sampling,  the  use  of  quasi- 
range  in  estimation  problems  is  satisfactory. 

Quasi-range  may  be  defined  as  follows.  Consider  the  ordered 

sample  values  Y, ,  Y.,  ...,  Y  where  a  <  Y  <  Y  <  ...<  Y  <  b.  If  the 
1  2  n  12  n 

r  smallest  and  r  largest  of  these  values  are  deleted,  the  range  of  the 
remaining  (n-2r)  values  is  defined  to  be  the  rth  quasi-range,  W^. 
Symbolically, 

W»Y  -  Y  . ,  n>  2r  +  1,  b-a>W  >0.  (1.1) 

r  n-r  r+1  r 

Note  that  sample  range  is  simply  the  quasi-range,  Wq. 

Many  authors  have  found  quasi-ranges  to  be  of  particular  interest 
and  have  developed  numerous  applications  for  them.  Among  others,  Harter 
[1959]  ,  Cadwell  [l953],  Chu  [1957]  ,  Benson  [1949]  and  Ghosal  [1957] 
have  derived  estimators  based  on  quasi-range,  and  Rider  [1959]  has 
studied  various  quasi-range  distribution. 

As  a  natural  extension  of  research  first  published  by  Gupta  [i960] 
regarding  the  distribution  of  order  statistics  from  a  gamma  density,  this 
study  of  quasi-range  is  intended  to  provide  further  information  needed  by 
the  analyst  in  applying  the  methods  of  simple  estimation  and  inference 
that  have  been  developed.  It  should  also  be  noted  that  Prescott  [1974 ] 
has  given  variances  and  covariances  of  order  statistics  from  the  gamma 
distribution.  Though  Tables  I,  II  and  III  show  the  moments  and  quantiles 
of  quasi-range  for  limited  sample  sizes  and  parameter  values  (n»l(l)  10  and 


a  •  0(1)  2)  only,  methods  described  here  can  be  easily  extended. 
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Distribution  of  Quasi-Range  in  Samples 


From  A  Gamma  Density 


The  density  function  of  the  standard  gamma  distribution  with 
parameter  a  is 


f(y)  -  (a!)"1  e"y  ya  (y  >  0) , 


(2.1) 


where  a  is  a  non-negative  integer.  The  corresponding  cumulative  distri¬ 
bution  function  may  be  expressed  as  a  partial  sum  of  Poisson  probabil¬ 


ities  : 


o  -y  i 

F(y)  -  1  -  Z  Z-JL  (y  >  0) . 
i*0  ' 


(2.2) 


Let  Wr  be  the  rth  quasi-range  from  a  sample  of  size  n  with  distri¬ 
bution  F(y) .  It  is  well  known  (Harter  [1959])  that  the  distribution  of 


U  can  be  expressed  in  integral  form  as 

r  b 


*  r!  (n-2r-2) Trl  /  “  F(y)1 


i-2r-2 


•  [1-F(y+Wf)]rf (y)f (y+w^)  dy  for  0  <  w^  <  b  -  a. 


(2.3) 


Substitute  (2.2)  into  (2.3)  and  let  a  (a,r)  be  the  coefficient  of  t 

-  ®  i  • 

in  the  expansion  of  (^E^tr/jJ]  (See  Prescott  [1974]),  then  (2.3)  can  be 
written  as 

r  i.  o*ii 

*(w_)  -  - 5 - “ - 2  1  (I  )<-*>  Z  ai 

(r!;  (n-2r-2)  I  (al;  i  -0  1  ij-0  l2  1 


n-2r-2 


z  ^n-2r-2j(_u  3  "z  3  (a>1  j  ?  (  4) 

.A  l  —A  J  l  —A  AC 


*3  alh 


i4 

r  4> 


i.-O  4 
4 


i5-0  5 


a(n-2r-2-i.)  a*r  i_  . 

Z  a  (a,n-2r-2-i.)  Z  a.  (a,r)  Z'(,7) 

i,"0  i6  3  i  -0  h  V 


i  (?,)  - 

l9*°  a 


(l4+i7«-l5-la-l9).e-“r'13+r+1) 


^  ~  (ij+n-r)  ‘  J 

•  (i2+i5+i6+ig+i9+a)l, w^  >  0. 


(i2+i5+i6+ig+i9+a+l) 


(2.4) 
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Ia  order  thac  this  expression  for  the  density  of  be  of  use,  it 

can  be  reduced  to  the  following  form, 

n-2r-2  -W  (r+j+1)  a(j+r+l)  . 

♦  (w  )  -  Z  e  Z  A(j,k)w  v  >  0,  (2.5) 

r  j-0  k-0  r  r 

where  A(j,k)  is  the  accumulation  of  all  coefficients  in  (2.4)  involving 

the  like  powers  of  wf  and  e.  This  function,  which  we  call  the  coefficient 

matrix  will  be  used  as  a  notational  convenience  throughout  this  paper. 


Although  it  is  impossible  to  express  A(j,k)  in  closed  form,  we  can  compute 
its  numerical  value  in  the  summation  process  of  (2,4)  by  evaluating  the 
coefficient  of  [axp{-w^ (r+j+1) }]  where  for  specified  values  of  n,  r  and  a 


k  ■  i.  +  i_  +  a 
4  7 


Using  ACl.k).  the  moment  generating  function  M(t)  associated  with 


quasi-range  can  be  written  as 


n-2r-2  a(j+r+l)  /,  w  , 

M(t)  -  E  Z  VJ  ...  ,  t  <  r  +  a  +  1  , 

j-0  k-0  (r+j+l-t)K  1 


(2.6) 


from  which  it  is  easily  seen  that  the  ith  moment  about  the  origin  of 
the  rth  quasi-range  is 


u|(n,r,a) 


n-2r-2  «(J+rfl)  A(1-k)(t 
j-0  k-o  (r+j+1) 


k+l+i  • 


(2.7) 


These  values  along  with  the  ith  central  moments  are  tabulated  in  Table  I. 

The  cumulative  distribution  function  for  can  be  obtained  from 
(2.5)  and  it  may  be  transformed  to  a  partial  sum  of  Poisson  prob¬ 


abilities  as 


n-2r-2  a(r+j+l)  k)k, 
*(w  )  -  Z  Z 

j-0  k-0  (r+j+1) 


M 

(  i-0 


k  [exp{-  w  (r+j+1))] (w  (r+j+1)} 

£  - - - 

(2.8) 


from  which  the  p-quantile  £  of  the  density  (2.5)  can  be  calculated  iter- 

P 

atively.  A  suitable  iterative  procedure  which  is  quite  efficient  in  com¬ 


puting  £  is  given  by 

p 
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(2.9) 


where  $(5p)  *  P  and  4> *  Pi  at  the  ith  Iteration.  To  initiate  the  pro¬ 
cess,  fairly  arbitrary  values  of  and  “a y  chosen  and  using  (2.8)  p^ 
and  p^  can  be  obtained.  These  values  are  tabulated  in  Table  III. 

Analysis  of  Results 

Ghosal  [1957]  has  published  results  of  an  analysis  of  quasi-range 
distributions  associated  with  samples  from  an  exponential  density. 

It  is  well  known  that  the  gamma  distribution  reduces  to  the  exponen¬ 
tial  distribution  when  a  *  0:  hence  a  limited  comparison  between 
Ghosal' s  data  and  that  of  this  study  is  possible.  Table  IV  presents 
this  comparison.  It  is  seen  that  u|(n,r,a)  and  u2(n,r,a)  from  Table  I 
agree  exactly  with  Ghosal 's  values,  but  values  from  Table  II  show  some 
disagreement.  Small  discrepancies  are  revealed  in  the  values  for-^,  and 
still  larger  errors  appear  for  To  find  the  cause  of  these  differ¬ 

ences,  Ghosal -s  formula  was  reapplied  with  more  accuracy  in  each  cal¬ 
culation.  Exact  agreement  between  data  sets  was  obtained  by  this  proce¬ 
dure,  leading  to  the  conclusion  that  Table  I  and  II  are  correct  for  o  *  0. 

Gupta  [1960]  has  published  tables  of  central  and  non-central  moments 
for  the  distribution  of  order  statistics  in  samples  from  a  gamma  density. 
Thus,  it  is  possible  to  form  a  direct  comparison  between  his  data  and 
that  of  this  study.  First  moments,  calculated  from  Gupta's  results, 
were  subtracted  from  those  of  Table  I  with  the  remainders  written  as 
absolute  deviations.  Results  of  this  comparison  are  presented  in  Table  V. 

A  maximum  deviation  between  the  two  data  sets  exists  when  n  ■  10,  a  •  2, 
and  r  «  3.  This  discrepancy,  amounting  to  0.00079,  reflects  a  deviation 
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from  Gupta's  results  of  0,06  percent  and  represents  the  largest  percent 
discrepancy  found.  Unfortunately,  no  recommendation  based  on  this  anal¬ 
ysis  can  be  made  regarding  the  relative  merit  of  either  data  set.  It 
can  only  be  surmised  that  discrepancies  may  be  due  solely  to 
machine  error  computation  (See  Prescott  [1974]).  Higher  moments  are  not 
compared  since  Gupta's  work  extends  only  to  the  derivation  of  formulas 
for  the  covariance  between  Y^  and  Y^.  Although  these  formulas  may  be  used 
to  calculate  quasi-range  variance,  the  amount  of  work  required  to 
obtain  a  comparison  is  considered  prohibitive.  To  obtain  still  higher 
moments,  original  deviations  for  covariance  between  powers  of  Yr  r  and 

Y. ,  are  required,  a  task  beyond  the  scope  of  this  study.  To  test  the 
r+i 

behavior  of  error  in  moments  of  higher  order,  an  alternative  computer  pro¬ 
gram  was  devised  that  differed  from  the  original  only  in  methods  used 
to  calculate  coefficients  A(j,k)  in  equation  (2.5).  Slightly  different 
values  for  P^(n,r,a)  were  computed,  but  in  general,  the  discrepancies 
tend  to  zero  as  the  order  of  the  moment  increased.  Two  examples  of 
this  phenomenon  appear  in  Table  VI. 

Another  independent  test  was  performed  to  determine  the  extent 
to  which  machine  error  accumulation  contributes  to  overall  error.  The 
basic  machine  program  used  for  all  computation  was  revised  so  as  to  per¬ 
form  calculations  in  double  precision  arithmetic;  only  small  differ¬ 
ences  less  than  0.0001  were  detected  while  computation  time  increased 
by  at  least  a  factor  of  two. 

In  addition  to  normal  machine  error,  values  in  Table  III  contain 
a  controlled  error  introduced  by  the  iteration  procedure  tolerance  limit. 
Thus  a  maximum  error  of  two  percent  can  exist  in  the  lower  tails  of  the 
distribution,  since  the  quoted  values  of  #(wf)  may  contain  as  much  as 
0.0001  error. 
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TABLE  I 


MOMENTS  OF  THE  QUASI-HAHGE  DENSITY 


a  =  0 


TABLE  I 


MOMENTS  OF  THE  QUASI-RANGE  DENSITY 


a  *  1 


n 

r 

u’ 

1 

u ' 

2 

u 1 

3 

u' 

4 

u 

2 

u 

3 

V 

4 

2 

0 

1.50000 

4.00000 

15.0000 

72.0000 

1.75000 

3.75000 

20.8125 

3 

0 

2.25000 

7.13889 

29.2500 

147.333 

2.07639 

3.84375 

24.0404 

4 

0 

2.74248 

9.70216 

42.3483 

221.445 

2.18098 

3.77775 

25.0139 

I 

0.77257 

1.09722 

2.21094 

5.70833 

0.50036 

0.59013 

1.73654 

5 

0 

3.10619 

11.6704 

54.3890 

293.197 

2.22195 

3.71378 

25.3323 

1 

1.28762 

2.34649 

5.48081 

15.5615 

0.68853 

0.66630 

2.42853 

6 

0 

3.39320 

13.7526 

65.5211 

362.346 

2.23880 

3.66255 

25.4053 

1 

1.67117 

3.57071 

9.26752 

20.3012 

0.77789 

0.70027 

2.7*526 

2 

0.52051 

0.50762 

0.70830 

1.27101 

0.23669 

0.19768 

0.40128 

7 

0 

3.62951 

15.4163 

75.8791 

428.936 

2.24497 

3.62173 

25.3746 

1 

1.97529 

4.7280b 

13.2999 

43.0566 

0.62631 

0.69622 

2.98760 

2 

0.91088 

1.18409 

1.97755 

4.01206 

0.35438 

0.25338 

0.63624 

a 

0 

3.8299b 

16.9145 

85.5740 

493.108 

2.24592 

3.58855 

25.2978 

1 

2.22641 

5.01167 

17.4336 

59.2309 

0.85479 

0.6.8824 

3.10655 

2 

1.22193 

1.91420 

3.63859 

8.11144 

0.42110 

0.27046 

0.76760 

3 

0.39246 

0.29211 

0.31307 

0.43306 

0.13807 

0.09004 

0.14036 

9 

0 

4.00373 

18.2741 

94.6958 

555.032 

2.2U20 

3.56093 

25.2013 

1 

2.43976 

6.82490 

21.5683 

76.4098 

0.67246 

0.67993 

3.18317 

2 

1. 47967 

2.65170 

5.56665 

13.3830 

0.46229 

0.27497 

0.66927 

3 

0.70647 

0.71717 

0.93851 

1.49405 

0.21806 

0.12374 

0.24225 

10 

0 

4.15673 

19.5212 

103.318 

614.375 

2.24107 

3.53740 

25.0974 

1 

2.62493 

7.77397 

25.7180 

94.2998 

0.88370 

0.67258 

3.22351 

2 

1.69910 

3.37648 

7.67506 

19.6406 

0.46953 

0.27455 

0.96067 

3 

0.96747 

1.20447 

1.82215 

3.23237 

0.26846 

0.13740 

0.31734 

4 

0.31500 

0.18973 

0.16538 

0.18663 

0.09050 

0.04860 

0.06167 
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TABLE  I 


MOMENTS  OF  THE  QUASI -RANGE  DENSITY 


1.67500  6.00000  26.2500  11*4.000 

2.61250  10.7676  52.0312  301.996 


3.1*251*1* 

0.97367 

3.67611 

1.62279 


0  1*.  23027 


2.10526 

0.65779 

1*.  52062 
2.1*6696 
1.15112 

1*. 76616 
2.80131* 
1.51*378 
0. 1*9666 

1*.  97860 
3.06781* 
1.66855 
0.891*11 

5.16550 

3.2961*7 

2.11*1*91* 

1.22357 

0.39892 


1U.7036 

1.69782 

18.0098 

3.61*635 

20.6718 

5.559U4 

0.79518 


75.9551 
1*.  11572 

96.0076 

10.3091* 

118.398 

17.5361* 

1.3536C 


U6O.20U 

12.571*9 

61U.61U 

31*. 8763 

763.993 
61*.  0997 
2.92363 


23.3970  137.351  906.102 


7.36737 

1.85765 

25.6583 

9.05773 

3.01072 

0.1*6066 

27.7073 
10.6355 
1*.  17317 
1.13315 

29.5816 

12.1102 

5.31528 

1.901*30 

0.3001*1* 


25.2509 

3.80536 

155.061 

33.1679 

7.03112 

0.60721* 

171.691* 

1*1.1229 

IO.7850 

1.82935 

167.384 

49.0226 

14.8951 

3.56196 

0.32382 


98.1885 

9.33924 

1047.05 

135.696 

19.0292 

1.02141 

1181.09 

175.614 

31.5621 


1310.56 

217.213 

46.4916 

7.73555 

0.44742 


TABLE  II 

KURTOSIS  AND  SKEWNESS  OF  THE  QUASI-RANGE  DENSITY 


2.00000 

6.00000 

1.60997 

4.08000 

1.46356 

3.48105 

2.00000 

6.00000 

1.33664 

3.19367 

1.49342 

3.44379 

1.33922 

3.02597 

1.28876 

2.63316 

2.00000 

6.00000 

1.30706 

2.91638 

1.17631 

2.24305 

1.45600 

3.23520 

1.28382 

2.83927 

1.10467 

2.01580 

1.22889 

2.34741 

2.00000 

6.00000 

1.26624 

2.76211 

1.05538 

1.86778 

1.10110 

1.91773 

1.43965 

3.14458 

1.25246 

2.73806 

1.01891 

1.76403 

1.01858 

1.66609 

1.20174 

2.21931 

2.00000 

6.00000 

3.T9592 
1.26467  I  2.57601 
1.17289 


1.75506 

1.05916 

0.83434 


1.21521 

1.05441 
0.8(5963 
0.80159  1.00888 
0.98776  1.1*0312 

1.781*79  1*. 52855 


2.93129 

1. 12746  1.94318 

5  1.71766 

66 


O.66165 

0.71036  0.76114 

0.90628  1.12565 

4.09185 
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TABLE  III 


QUANTILES  OF  THE  QUASI -RANGE  DENSITY 


0.005 

0.010 

0.025 

0.050 

0.00500 

0.01002 

• 

0. 02532 

0.05126 

0.07301 

0.10538 

0.17216 

0.25311 

0.16740 

0.24300 

0.34585 

0.459 69 

0.00250 

0.00501 

0.01266 

0.02564 

0.30863 

0.38007 

0.50671 

0.64030 

0.04213 

0.06072 

0.09908 

0.14544 

0.42479 

0.50748 

0.64990 

0.79688 

0.11740 

0.15207 

0.21580 

0.28596 

0.00167 

0.00334 

0.00844 

0.01709 

0.53375 

0.62354 

0.77613 

0.93379 

0.20479 

0.25108 

0.33334 

0.41944 

0.02979 

0.04291 

0.07000 

0.10271 

0.63195 

0.72900 

0.89250 

1.05498 

0.29324 

0.34643 

0,44390 

0.54159 

0.08629 

0.11183 

0.15655 

0.20991 

0.00125 

0.00251 

0.00633 

0.01282 

0.72437 

0.82553 

0.99556 

1.16351 

0.37824 

0.44069 

0.546y0 

0.65251 

0.15514 

0.18993 

0.25139 

0.31653 

0.02307 

0.03323 

0.05421 

0.07952 

0.81127 

0.91547 

1.08982 

1.26206 

0.45961 

0.52729 

0.64153 

0.75352 

0.22698 

0.26926 

0.34294 

0.41730 

0.06841 

0.08863 

0.12568 

0.16633 

0.00100 

0.00200 

0.00506 

0.01026 

0.100 


0.10529 

0.38013 

0.62398 

0.05264 


0.21791 

0.99684 

0.36639 

0.03510 

1.1434? 

0.53824 

0.15381 

1.27163 

0.67272 

0.28325 

0.02632 

1.36506 

0.79261 

0.40525 

0.11905 

1.48857 


0.51714 

0.22431 

0.03106 


0.500 


0.69307 
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TABLE  III 

QUANTILES  OF  THE  QUASI-RANGE  DENSITY 


q  »  0 


n 

11 

!  P 

B 

0.900 

0.950 

0.975 

0.990 

0.995 

2 

0 

2.30240 

! 

2.99567 

3.60498 

4.60437 

5.29340 

3 

0 

2.96973 

3.67609 

4.37217 

5.29507 

5.96538 

4 

0 

3.36629 

4.07731 

4.77571 

6.39047 

1 

1.15120 

1.40784 

1.84240 

2.64670 

5 

0 

3.61*978 

4.36286 

5.06246 

5.96700 

6.67798 

1 

1.63025 

1.99987 

2.35957 

2.83150 

3.16206 

6 

0 

3.86992 

4.56282 

5.26502 

6.20091 

6.90103 

1 

1.91*796 

2.32575 

2.69276 

3.16979 

3.52290 

2 

0.7671*7 

0.99856 

1.22832 

1.53478 

1.76445 

7 

D 

1*. 04991 

4.76439 

5.46700 

6.39209 

7.08331 

H 

2.16698 

2.57036 

2.94073 

3.42061 

3.77512 

2 

1.13774 

1.39188 

1.63814 

1.95970 

2.19767 

8 

4.20322 

4.91802 

5.62091 

6.54613 

7.23741 

n 

2.37942 

2.76620 

3.13630 

3.62038 

3.97576 

2 

1.39982 

1.66400 

1.91677 

2.24750 

2.46940 

3 

0.57560 

0.74892 

0.92126 

1.15113 

1.32342 

9 

0 

4.33661 

5.05119 

5.75431 

6.67967 

7.37111 

1 

2.54087 

2.92966 

3.30377 

3.78656 

4.14252 

2 

1.60483 

1.87563 

2.13497 

2.46779 

2.71177 

3 

0.87661 

1.07119 

1.25935 

1.50448 

1.68548 

0 

4,45434 

5.16867 

5.87193 

6.79731 

7.48870 

1 

2.67921 

3.07020 

3.44561 

3.92966 

4.20716 

2 

1.77394 

2.04899 

2.31117 

2.64607 

2.89041 

3 

1.09902 

1.30390 

1.50080 

1.75395 

1.93967 

U 

0.46045 

0.59908 

0.73688 

0.92057 

1.05808 
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QUANTILES  OF  THE  QUASI -RANGE  DENSITY 


a  =  1 


m 

P 

B 

0.005 

0.010 

0.025 

0.050 

0.100 

0.500 

■ 

0 

0.01003 

0.02000 

0.Q5000 

O.lQQQl* 

0.20122 

1.11*619 

9 

KB 

0.15027 

0.21312 

0.33959 

0.1*8555 

0.70285 

1.96018 

D 

0.36050 

0.1*81*75 

0.66625 

0.85909 

1.12332 

2.1*7699 

O.OOULo 

0.00896 

0.02256 

O.OU55I* 

0.0928U 

0.5731*6 

5 

§:'U 

0.61465 

0.71*251* 

0. 95720 

1.17211 

1.1*5780 

2.85266 

0.07813 

0.11172 

0.18037 

0.2611*9 

0.381*91 

1.12388 

6 

x;.:  ■ 

0.63269 

0.97362 

1.20669 

1.1*31*11* 

1.73171 

3.11*683 

1 

0.22091 

0.28193 

0.39321 

0.51212 

0.6771*3 

1.52622 

2 

0.00287 

0.00575 

O.OIU50 

0.02932 

0.05996 

0.37976 

7 

0 

1.031*10 

1.17771* 

1.1*2216 

1.65722 

1.96137 

3.38761 

1 

0.38293 

0.E6LL2 

0.601*13 

0.7U625 

0.93501* 

1.81*107 

2 

0.05267 

0.07579 

0.12231* 

0.17088 

0.261*88 

0.79123 

0 

0 

1.20573 

1.36039 

1.61051+ 

1.6501*7 

2.15960 

3.59110 

1 

0.5^697 

0.61*110 

0.80031* 

0.95822 

1.15963 

2.09921* 

2 

0.15572 

0.19973 

0.28006 

0.36673 

0.1*871*0 

1.11605 

3 

0.00211 

0.001*23 

0.01062 

0.02160 

0.01*1*23 

0.28331* 

9 

D 

1.36572 

1.52221 

1.77725 

2.0201*3 

2.33252 

3.76710 

m 

0.70230 

0.80569 

0.97215 

1.11*357 

1.35663 

2.31730 

2 

0.26009 

0.31*052 

0.1*1*512 

0.55203 

0.69U79 

1.38255 

3 

0.03998 

0.05737 

0.09316 

0.13591* 

0.20188 

0.61061* 

10 

Q 

1.U9910 

1.66722 

1.92535 

2.17261 

2.1*8593 

3.92191* 

0.81*761* 

0.95677 

1.13903 

1.31120 

1.53102 

2.50571 

H 

0.1*1212 

0.1*8390 

0.60595 

0.72609 

0.88332 

I.608U9 

H 

0.12051* 

0.151*96 

0.21776 

0.23582 

0.38106 

0.68207 

1 

0.00167 

0.00335 

O.OO8U5 

l— — 

0.01710 

0.03503 

0.22581 
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TABLE  III 


QUANTILES  OF  THE  QUASI -RANGE  DENSITY 
_  a  a  1 


P 


m 

0.900 

0.950 

0.975 

0.990 

0.995 

0 

3.27175 

4.11295 

4.92794 

5.90947 

6.77216 

0 

U. 16697 

5.00322 

5.31462 

6. 86676 

7.64312 

0 

4.69067 

5.52970 

6.33657 

7.36306 

6.15586 

1 

1.72306 

2.17995 

2.62111 

3.19314 

3.61311 

0 

5.07275 

5.90270 

6.70622 

7.74376 

8.51912 

1 

2.40290 

2.67273 

3.32165 

3.69623 

4.31575 

0 

5.30439 

6.19068 

6.99130 

6.03133 

6.79966 

1 

2.6^012 

3.32395 

3.77413 

4.34775 

4.76602 

2 

1.17316 

1.49210 

1.80047 

2.20036 

2.49371 

0 

5.60113 

6.42526 

7.22407 

8.26119 

9.02326 

1 

3.10566 

3.66056 

4.11055 

!  4.60296 

5.10014 

2 

1.71395 

2.05146 

2.37266 

2.78129 

3.07833 

5.60066 

6.62250 

7.41953 

6.45465 

9.22049 

B 

3.45^06 

3.92910 

4.37847 

4.94955 

5.36566 

2 

2.09308 

2.43834 

2.70309 

3.1753d 

3.47334 

3 

O.60966 

1.13567 

1.37447 

1.66435 

1.91177 

0 

5.97160 

i  6.79259 

7.56613 

6.62159 

9 . 36646 

1 

3.67748 

4.15241 

4.60109 

5.17104 

5.5o644 

2 

2.38878 

2.73751 

3.06497 

3.47729 

3.77511 

3 

1.33711 

1.60321 

I.650IO 

2.17697 

2.40930 

0 

b. 12258 

6.94204 

7.73634 

6. 76o52 

9.53261 

1 

3.86839 

4.34273 

4.79024 

5.35765 

5.77002 

2 

2.63223 

2.96334 

3.31366 

3.73052 

4.03501 

3 

1.66430 

1.93972 

2.19844 

2.52363 

2.75732 

4 

0.71653 

0.91741 

1.11255 

1.36651 

1.55314 

k 


TABLE  III 

QUANTILES  OF  THE  QUASI -RANGE  DENSITY 


a*2 


■ 

p 

B 

0.005 

0.010 

0.025 

0.050 

0.100 

0.500 

e 

0 

0.01336 

0.02667 

0.06667 

0.13337 

0.26772 

1.48193 

H 

0 

0.20333 

0.26564 

0.45424 

0.04757 

0.93238 

2.50734 

D 

0 

0.50879 

0.61*998 

0.89053 

1.14226 

1.48274 

3.15140 

1 

0.00585 

0.01170 

0.02945 

0.05941 

0.12098 

0.73783 

5 

0 

0.82902 

0.99379 

1.27530 

1.55330 

1.91741 

3.61662 

1 

0.10271* 

0.11*676 

0.23661 

0.34240 

0.50252 

1.43883 

6 

0 

1.11818 

1.300U1 

1.60341 

1.89431 

2.27022 

3.97875 

1 

0.28977 

0.37138 

0.51652 

0.67068 

0.88390 

1.94767 

2 

0.00372 

0.0071*5 

0.01877 

0.03794 

0.07754 

0.48703 

7 

0 

1.37959 

1.56939 

1.88513 

2.13369 

2.56482 

4.27413 

1 

0.50378 

0.61187 

0.79431 

0.97691 

1.21834 

2.34410 

2 

0.06891 

0.09870 

0.15980 

0.23238 

0.34338 

1.01193 

6 

D 

0.02391 

0. 01+771 

0.11927 

0.23854 

0.47830 

3.19375 

m 

0.72115 

0.81*439 

1.05052 

1.25044 

1.50886 

2.66761 

2 

0.20387 

0.26082 

0.36499 

0.47662 

0.63212 

1. 42440 

3 

0.00272 

0.00546 

0.01376 

0.02784 

0.05696 

0.36266 

9 

D 

0.02525 

0.05044 

0.12609 

0.25216 

0.50493 

3.37813 

m 

0.01676 

0.03332 

0.08329 

0.16657 

0.33603 

2.42899 

2 

0.36761 

0.44620 

0.58122 

0.71376 

0.90124 

1.76219 

3 

0.05181 

0.07431 

0.12057 

0.17576 

0.26060 

0.78040 

10 

EE 

0.02654 

0.05303 

0.13262 

0.26519 

0.53041 

3.54876 

B 

0.01752 

0.03496 

0.08741 

0.17478 

0.35075 

2.59041 

n 

0.011*26 

0.02803 

0.07006 

0.14066 

0.29397 

1.94924 

H 

0.15658 

0.20113 

0.28240 

0.37020 

0.49250 

1.12607 

1 

0.00215 

0.00430 

0.01086 

0.02197 

0.04500 

0.28867 

TABLE  III 


} 


QUANTILES  OF  THE  QUASI -RANGE  DENSITY 


o  =  2 


B 

p 

0.900 

0.950 

0.975 

0.990 

0.995 

H 

9 

U. 01037 

1*.  96355 

5.03349 

7.05906 

7.91742 

9 

0 

5.07972 

6.02017 

6.92124 

a. 07669 

0.92230 

9 

0 

5.71162 

6.61*201 

7.53313 

0.67744 

9.51617 

1 

2.15116 

2.69395 

3.21071 

3.07134 

4.35110 

5 

0 

6.15963 

7.06213 

7.96644 

9.10309 

9.93710 

D 

2.9*3259 

3.53200 

4.05055 

4.70579 

5.17986 

6 

M 

6.50574 

7.1*2170 

6.30004 

9.43101 

10.2625 

3.53763 

4.07710 

4.50265 

5.24297 

5.71319 

1.1*7396 

1.66055 

2.22963 

2.70262 

3.04506 

T 

0 

6.76563 

7.69744 

0.57246 

9.69697 

10.5269 

D 

3.93606 

4.43300 

4.99551 

5.64059 

6.10667 

9 

2.11*266 

2.54566 

2.92510 

3-40241 

3.74659 

fi 

R 

1*. 79530 

5.01436 

5.13140 

5.20337 

5.22781 

9 

4.26202 

4.00662 

5.31629 

5.95746 

6.42078 

2 

2.61016 

3.01350 

3.39975 

3.37602 

4.21904 

3 

1.12121 

1.42296 

1.71239 

2.05399 

2.35403 

9 

D 

u. 97991* 

5.19647 

5.31130 

5.36197 

5.40594 

n 

3.60b93 

3.00306 

3.91344 

3.90470 

4.00965 

2 

2.971*71 

3.38592 

3.76904 

4.25049 

4.60129 

3 

1.67762 

1.99091 

2.30077 

2.67040 

2.94761 

10 

0 

5.11*91*7 

5.36306 

5.47697 

5.54731 

5.57095 

BE 

3.71*002 

3.92059 

4.02941 

4.09351 

4.11570 

3.06666 

3.30293 

3.59500 

3.76396 

3.83304 

9 

2.06667 

2.41903 

2.73064 

3.12633 

3.42082 

1 

O.9OU02 

1.15297 

1.39213 

1.7009b 

1.92703 

TABLE  V 


COMPARISONS  WITH  GUPTA'S  DATA 
Absolute  Deviations 


n 

r 

a  *  0 

a  =  1 

e  =  2 

2 

0 

- 

- 

- 

3 

0 

- 

- 

- 

4 

0 

- 

0.00001 

0.00001 

1 

- 

- 

- 

5 

0 

- 

- 

- 

I 

- 

0.00001 

- 

6 

H 

- 

- 

- 

n 

mm 

0.00001 

0.00001 

- 

0.00001 

7 

H 

mm 

• 

0.00001 

• 

- 

0.00001 

6 

0.00001 

0.00002 

n 

mm 

- 

0.00003 

2 

0.00002 

0.00002 

3 

- 

0.00001 

0. 00003 

9 

p 

- 

0.00001 

0.0G002 

fl 

- 

0.00002 

0.00003 

2 

- 

0.00004 

0.00006 

3 

- 

0.00001 

0.00004 

10 

B 

- 

0.00001 

0.00004 

0.00001 

0.00003 

0.00018 

0.00002 

0.00001 

0.00023 

B 

- 

0.00007 

0.00079 

1 

- 

0.00012 

0.00028 
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TESTS  OF  ACCURACY  IW  HIGHER  MOMENTS 


yj(n,r,a) 


D 

Primary 

Program 

Alternate 

Program 

1  ■ 

D 

1.69910 

1.69922 

3.37646 

3.37656 

B 

7.67506 

7.67516 

HU 

H 

19.6406 

19.6407 

B 

2.14494 

2.14516 

o 

rA 

II 

a 

5.31526 

5.31542 

<x  =  2 

14.8951 

14.8952 

r  =  2 

B 

46.4916 

46.4915 
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Appendix  B 


ON  ESTIMATING  THE  SCALE  PARAMETER  OF  THE 
RAYLEIGH  DISTRIBUTION  FROM  DOUBLY  CENSORED  SAMPLES 


Introduction 


This  paper  is  concerned  with  estimating  the  scale  parameter  of  the  Rayleigh 
distribution  from  censored  samples.  The  Rayleigh  distribution  arises  as  a  con¬ 
sequence  of  finding  the  resultant  amplitude  of  several  coplanar  random  amplitude 
vectors  which  are  normally  distributed  (Siddiqui  [1962]).  Therefore  it  is  useful 
in  the  analysis  of  acoustic  data  or  other  data  obtained  from  measurements  of 
amplitudes  of  electromagnetic  waves  received  through  a  scattering  medium. 

This  distribution  is  also  useful  in  communication  engineering.  Since  for 
many  reasons  the  samples  could  be  censored,  it  is  important  to  consider  the 
analysis  of  such  data. 

The  Rayleigh  distribution  is  characterized  by  the  probability  density  func¬ 
tion  (p.d.f.) 

((2x/K)exp{-x  /K) ,  0  <  x  <  “ 

fW  *  {  (l.l) 

(  0  ,  elsewhere 

for  positive  values  of  K,  with  expectation  /KtF/2  and  variance  K(1-tt/4).  To 
estimate  the  parameter  K,  we  employ  the  methods  used  by  Tiku  [1967,  1967(a), 
1968,  1968(a),  1968(b)]  for  the  censored  samples  from  normal,  exponential, 
logistic,  log-normal  distributions  and  progressively  censored  samples  from 
normal  distribution  respectively,  as  described  below. 

The  p.d.f.  given  in  (1.1)  can  be  reduced  to 

(2z  exp(-z2)  ,  0  <  z  <  - 

(1.2) 

0  ,  elsewhere, 

where  x  *  Z//K 

Let  g(z)  -  fz(z)  /Fz(z),  where  Fz(z)  is  the  probability  integral  of  fz(z) 
given  in  (1.2).  Then,  over  a  small  interval,  say  a  <  x  <  b,  the  linear 
approximation  a+  B  z  to  g(z)  is  a  reasonable  one,  where  a  and  B  are  constants 


-- 


8  -  (g(b)  -  g(a) }/ (b-a) 
a  -  g(a)  -  a8  . 

Some  numerical  comparison  between  g(z)  and  a  +  8z  for  various  sample  sizes 
are  given  in  the  example.  The  substitution  a  +  8z  for  g(z)  in  the  likelihood 
equation  results  in  a  solution  (K^)  which  is  easy  to  compute  and  which  is 
asymptotically  equivalent  to  the  maximum  likelihood  estimator  (MLE)  accord¬ 
ing  to  Tiku's  results  [1967,  1967(a),  1968,  1968(a),  1968(b)],  To  obtain 
greater  accuracy,  i.e.,  to  get  an  approximation  which  is  practically  same 
as  an  actual  MLE,  it  is  suggested  to  try  a  linear  approximation  once  again 
instead  of  using  expensive  and  time  consuming  iteration  procedures. 

In  the  remainder  of  this  paper  the  maximum  likelihood  equation  for 
finding  K  will  be  set  up  and  solved  (for  a  doubley  censored  sample), 
expressions  for  the  bias  and  variance  of  Kc  will  be  developed,  a  numerical 

example  will  be  given  and  an  improved  estimator  of  K  (namely  K')  will  be  made 

c  ^ 

Derivation  of  the  Estimator  K 

c 

Let  Xj ,  X^,  ...,  Xn  be  a  random  sample  from  the  Rayleigh  distribution 
with  the  smallest  ®  n  and  the  largest  r 2  «  q2n  observations  being  cen¬ 
sored,  where  q^  and  q2  are  fixed  and  decided  in  advance  (Type  II  censoring). 
The  remaining  sample  values,  arranged  in  order  of  magnitude,  are 

Y  Y  Y  Y 

rj+1  rj+2*  n-r2-l*  n-r2’ 

i.e.  forming  a  doubly  censored  Rayleigh  sample  of  size  n  "  "  r2*  T^ie 

probability  density  function  of  this  censored  sample  (see.  for  example 


Saw  [1961])  is 


— 7(2//K)a  ^  r2  fY2  z  jexp  j  -l  2  zZ  \ 
rlJr2f  Li-r^y  (  i-r^l  M 

where  z^  ■  y^//K. 


Ks,+1)  fl 


(2.1) 
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Taking  logarithms  of  (2.1)  and  denoting  it  by  L  yields: 


logfr-TTT  1  *  Cn-r.-r,)  log  (2//K)  ♦  l  log  z 
(rl ‘ r2  *  f  1  i=r1+l  1 


n-r„ 


(2.2) 


1  zT  +  r.  {log  F_ ( 2  .)}  -  r.  i 

.,i  1  Z  rl+1  2  n_r2 


i=r^+l 


Taking  the  partial  derivative  of  L  with  respect  to  K  and  simplifying 


3L 

5k 


2 

_  _  n-r_  Pi,Z  . f _ f 2  .)  r_z 

,  rw"  i  t  2  z2  \ 1  tj1]  i  Txi  ,  2n~r^  f2 

L  K  J  K  l  i  L  2K  J^C,  )  +  K  -t2- 


3) 


Setting^  equal  to  zero  in  equation  (2.3)  and  solving  for  K  would  give  the 
ordinary  maximum  likelihood  estimator.  However,  this  is  a  complicated  nonlinear 
equation  because  of  the  term 


g  (2r  +1)  =  f2(zr  +1)/Fz(zr  +l) 


2zr1+l  exp  <"*r1+l> 


1  *  1-  exp  (-z2i+1) 


(2.4) 


which  is  implicit  in  K,  thus  precluding  any  exact  solution  to  (2.3).  Hence  an 
iteration  procedure  seems. to  be  the  only  solution  to  this  problem,  but  it  is 
not  only  time-consuming  and  expensive  for  computation  but  also  sometimes  diffi¬ 
cult  to  converge  as  is  the  case  with  the  Rayleigh  distribution*.  Hence  instead 
of  an  iteration  procedure,  we  are  proposing  that  g(zf  +1)  be  replaced  by  a  linear 
approximation 
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V1 


)  =  a  S  z 


V1' 


(2.5) 


Consider  now 


*We  experienced  divergences  even  with  the  actual  population  parameter  as  an  ini 
tial  guess  when  censoring  is  a  little  bit  heavy. 
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3L  _  3 L*  _  rl+r2-n 
3K  ‘  3K  K 


+  K 


n-r„ 


i=rl+1 


’Vr  ♦ll  r  1  r2zn-r. 

-  hH  1°  * fcv  J  * 


(2.6) 


In  this  equation  a  and  8  are  such  that 

8  =  (gz(h2)  -  gz(h1)>  /Ch2-h1)  and 


a  =  Sz  Chi )  -  hj  8 


(2.7) 


where  g  is  given  by  (2.4)  and  the  interval  (h^.h^)  is  chosen  in  such  a  way  that 
z  .  is  sufficiently  close  to  h,  or  h_.  But  a  difficulty  arises,  because  we 

Tj+1  1  2 

don't  know  the  exact  point  of  z^  since  it  depends  on  the  parameter  K  we  are 
going  to  estimate.  This  difficulty  can  be  eliminated,  for  sufficiently  large 
n-r^-r2,  by  choosing  h^  and  h^  so  that 


W  - q! 


q1(i-q1) 


(2.8) 


w  ■ 


The  reasoning  behind  this  choice  of  interval  end  points  is  that  it  is  logical 
to  think  of  z^  ^  as  an  estimate  of  the  point  below  which  100*q^  percents  of 
the  population  represented  by  fz(;)  lies.  Also,  the  secant  between  h^  and  h2  has 
smaller  maximum  error  than  the  tangent  at  h  or  than  the  secant  through  hj'  and 
h  '  if  one  recalls  that  the  probability  is  small  that  the  Z  .  will  fall  out- 
side  the  interval  h  '  and  h2',  where  Fz(h)  =  q^, 


Fz(h1')  *  qj  -  3^-q^l-qp.  ,  (2.9) 

and  Fz(h2’)  ■  ♦  3Vn  * 


(See  figure  1) 
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From  (2.8)  it  can  be  seen  that' as  n  become  large,  F^Ch^)  approaches 

from  the  left  at  the  same  rate  that  F_(h'j  approaches  from  the  right.  Thus 

the  interval  (h^,h2)  shrinks  to  a  single  point,  and  a  and  3  can  be  obtained 
by  evaluating  the  derivative  of  gz(z)  at  the  point h,  and  in  this  case  all  3 
lines  in  the  figure  coincide.  The  degree  of  accuracy  of  the  linear  approxi¬ 
mation  is  related  to  the  width  of  the  interval:  The  smaller  the  width  of 

(hj,h2),  the  smaller  the  error  of  the  approximation  which  is  obvious  from  the 

figure.  But  one  disadvantage  is  that  the  decrease  of  interval  width  f^-h^  is 
rather  slow  as  n  increases  though  the  approximation  war  reasonably  good  when 
n  is  as  small  as  10  (see  §4) .  Hence  it  is  suggested  in  55  that  a  possible 
acceleration  of  the  approximation  be  used. 

3L  * 

Equation  (2.6)  can  now  be  solved  analytically  by  setting  s  0  and 

oK 

carrying  out  the  algebra.  After  substitution  of  y  *  /ff  and  some  simplifi¬ 
cation,  it  follows  that 


-—•(GK+B-D/K  )  -  0, 
21T 

where  G  •  2(r^+r2~n) 


B 


n-r2  , 

2  t  yt  +  2 

i-r,+l 


r2y 


2 

n-r2 


rlyr1+ie 


D  ■  r.y  ,a. 
I'r.+l 


(2.10) 
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In  order  to  solve  equation  (2.10)  for  K  >  0, 

set  T  ■  and  rewrite  it  as  a  quadratic  equation  in  T  as 

GT2  -  DT  +  B  »  0.  (2.H) 

Since  4c  is  positive  we  take  the  positive  root  of  the  equation  (2. 11 )  as 

T,  i.e.  _ _ 

D  -  -4GB 
T - 2G -  * 

Here  note  that  G  is  negative  and  B  and  D  are  positive.  Hence  the  estimator 
K  is  the  square  of  T  as 

Kc  *  |  (°2"2GB^  "  %^-4Gb|/2C2  (2.12) 

Properties  of  the  Estimator  Kc 

The  estimator  K£  is  same  as  an  MLE  except  that  a  linear  approximation 
was  used  to  solve  the  likelihood  equation.  It  is  expected  that  the  properties 
of  Kc  should  be  similar  to  those  of  an  MLE.  While  calculation  of  the  expected 
value  of  K£  from  equation  (2.12)  would  be  extremely  difficult,  we  can  discuss 
the  approximate  conditional  bias  of  K£  in  the  asymtotic  case.  Following  Tiku 
[1967,  p.  160],  we  have  the  approximate  conditional  bias  given  by 
B1  3  EC^’)/r2CK) 

where  (3.1) 

2  32L* 

R  (K)  *  -EC—,) 

3  K" 

for  large  values  of  n  -  r^  -  Tj. 

The  bias  can  be  calculated  from  the  following  equation: 

n-r. 

,  rr.+r^-ni  2  - 

E(||r  >  “  I  1  *  (1/K)  I  E(zf )  -  (r  a/2K)  E(Z  ) 

.  i-r^l  1  1  V1 

-  (ri0/2K)E(zJ  )  ♦  (r2/K)  EfZ2^).  (3.2) 

Expected  values  of  the  order  statistics  from  the  Rayleigh  distribution  are 
given  by  (see  Sarhan  and  Greenberg  [1962]) 


and 


E(Z.) 


(i-l)  Un-i) !  j0  [— ^p7?J 


(3.3) 


Differentiating  (2.6)  with  respect  to  K  and  taking  expectations  and  inserting 
the  minus  sign  gives 


R2(K) 


-E  (&  ) 

\  n/ 

T~  J  E<ZP - T 

K  i=r, +1  1  4K  1 


n-r. 

r,*r2-n  „  2 


3r,ct 


rie 


2r. 


-  "T1  E(Zr  ♦  !>  *  — T 
K  1  K 


f2-J 


(3.4) 


where  expectations  are  also  given  by  (3.3). 

Also  the  asymptotic  variance  of  can  be  obtained  from  (3.4)  by  use  of  the 

asymptotic  property  (See  Kendall  and  Stuart  [1961]): 

Var(K)  =  -rr—  (3.5) 

r(K) 


As  Tiku  [1967]  justified,  since  Kc  is  an  asymptotic  MLE,  its  asymptotic  properties, 
in  no  doubt,  are  the  same  as  MLE.  We  can  show  that  numerically,  though  they  are 
complex,  the  bias  in  (3.1)  would  be  zero  when  n  is  large.  Furthermore,  the 
attractiveness  of  the  estimator  is  enhanced  by  the  fact  that  it  can  be  computed 
easily  without  worrying  about  divergence,  without  having  to  resort  to  iterations 
or  procedures  requiring  expectations  of  order  statistics  (see  Lloyd  [1952]). 

Table  1  below  shows  the  values  of  and  VarfK^)  for  various  censoring  when  the 
sample  size  is  10  and  K  ■  3. 
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Table  1.  Asymptotic  bias  and  variance  of  the  estimator  Kc  for  n  ■  10  with 
proportion  censored  from  below  and  ^2  abovo • 


*1 

*2 

Cl) 

(2) 

.1 

.2 

0.019 

0.469 

.1 

.4 

-  .092 

.160 

.1 

.6 

-  .087 

.203 

.1 

.7 

-  .OSS 

.262 

.2 

.2 

-  .009 

.687 

.2 

.4 

-  .167 

.265 

.2 

.6 

-  .173 

.363 

.3 

.4 

-  .217 

.279 

.3 

.5 

-  .242 

.403 

.4 

.4 

0 

-  .055 

.170 

(1)  =  Bias/o 

(2)  =  Var(Kj/o2, 

where 

a2  =  Variance  of  the  Rayleigh  distr.  =  K[l-ir/4) 

Numerical  Example 

Using  .the  closed-form  cumulative  distribution  function  F(x)  for  the  Rayleigh 
distribution  and  the  probability  integral  transformation,  a  random  sample  of  100 

observations  was  generated  from  the  Rayleigh  distribution  with  the  population 
parameter  K  ■  3.000.  The  sample  of  size  10  and  30  were  chosen  from  the  original 

•ample  of  size  100. 

The  following  analysis  of  data  is  given. 

4.1  Linear  approximation  a  +0  zr  to  +J)  . 

The  following  table  shows  the  linear  approximation  a+0  z  to  g(z  )  for 

r^  +  1  tj  ^ 1 

various  sample  sizes  with  proportions  qj  and  q,  censored  from  the  left  and  the 
right,  respectively. 
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Table  2.  n  *  10 


ql 

q2 

hl 

h2 

a  ♦  eCy^/'T) 

.1 

.1 

.071 

.465 

7.09 

IS. 65 

.1 

.2 

6.33 

13.82 

.1 

.3 

6.33 

13.82 

.2 

.1 

.276 

.628 

4.16 

4.98 

.2 

.2 

3.67 

4.39 

.2 

.3 

3.67 

4.39 

.3 

.1 

.410 

.767 

4.16 

4.28 

.3 

.2 

3.67 

3.94 

.3 

.3 

3.67 

3.94 

.4 

.1 

.530 

.899 

V 

3.66 

3.49 

.4 

.2 

3.19 

3.22 

Ke 

is  an  actual 

.  MLE  by 

iteration 

• 

i  3. 

n  =  30 

qi 

q2 

hl 

h2 

sfVi/y^ 

*  * 

.1 

-.1 

.215 

.410 

8.65 

8.84 

.1 

.2 

7.78 

8.27 

.1 

.3 

7.20 

7.83 

.2 

.1 

.368 

.564 

4.S6 

4.69 

.2 

.  .2 

4.04 

4.24 

.2 

.3 

3.71 

3.90 

.3 

.1 

.494 

.695 

4.06 

3.89 

.3 

.2 

3.58 

3.58 

.3 

.3 

3.27 

3.34 

.4 

.1 

.609 

.820 

2.83 

2.79 
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Table  4.  n  ■  100 


*1 

q2 

hl 

h2 

g(z) 

a  ♦  8  z 

.1 

.1 

.269 

.373 

7.52 

7.41 

.1 

.2 

7.10 

7.11 

.1 

.3 

7.22 

7.20 

.2 

.1 

.417 

.523 

4.13 

4.17 

.2 

.2 

3.87 

3.92 

.2 

.3 

3.94 

3.99 

.3 

.1 

.541 

.651 

2.95 

2.97 

.3 

.2 

2.74 

2.76 

.3 

.3 

2.80 

2.83 

.4 

.1 

.657 

.772 

2.20 

2.21 

.4 

.2 

2.02 

2.03 

It  can  be  seen  from  our  earlier  discussion,  the  approximation  is  better  when 
n  is  large  or  hj-h^  is  small  but  when  the  sample  size  is  as  small  as  10,  the 
approximation  is  not  good;  hence  an  improvement  is  suggested  in  55. 


4.2  Comparison  With  an  actual  MLE  by  iteration. 

Table  5  shows  comparison  between  and  an  MLE  by  iterations*  for  n  a  30. 
Table  5. 


ql 

q2 

K 

c 

MLE 

Error 

.1 

.1 

2.934 

2.941 

.002 

.1 

.2 

2.380 

2.399 

.007 

.1 

.3 

2.0S4 

2.078 

.011 

.1 

.4 

1.887 

1.917 

.015 

.1 

.5 

1.915 

1.951 

.018 

.2 

.1 

2.938 

2.956 

.011 

.2 

.1 

2.390 

2.417 

.012 

.2 

.3 

2.071 

2.097 

.013 

.2 

.4 

1.914 

1.940 

.016 

.2 

.5 

1.946 

1.978 

.015 

.3 

.1 

2.985 

2.951 

.011 

.3 

.2 

2.411 

2.410 

.000 

.3 

.3 

2.076 

2.090 

.006 

.3 

.4 

1.909 

1.931 

.011 

.3 

.5 

1.942 

1.967 

.012 

.4 

.1 

3.003 

2.990 

.004 

.4 

.2 

2.432 

2.452 

.007 

.4 

• 

.3 

2.10S 

2.132 

.012 

.4 

.4 

1.951 

1.981 

.015 

.4 

.5 

1.991 

2.029 

.018 

Error  *  |K  -MLE) /MLE 

V 


'Iterated  approximation  is  obtained  by  Newton- Raphson’s  method  (See,  for 
example,  Dahlquist  [1]). 
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Improvement  by  linear  Approximation  Twice 
Sometimes  a  greater  accuracy  for  an  approximation  of  an  MLE  is  needed  and 
as  it  was  seen  in  §4,  our  estimator  K  is  not  sufficient  for  that  purpose 
when  the  sample  size  is  smaller  than  30.  To  achieve  greater  accuracy,  one 
solution  is  to  use  Kc  as  an  initial  guess  and  try  an  iteration  procedure, 
which  will  give  faster  convergence  in  every  case.  But  considering  the  compu¬ 
tation  time  involved  in  iteration,  it  is  suggested  to  try  one  more  linear 
approximation  using  Kc  as  an  estimator  of  K,  which  resulted  in  an  estimator 
Kc'  and  showed  almost  the  same  accuracy  (agreed  with  an  actual  MLE  by  itera¬ 
tion  to  two  decimal  digits  when  n  *  10)  as  an  MLE  by  iteration.  The  second 


linear  approximation  procedure  is  as  follows:  Since  an  estimator  Kc  of  an 

actual  MLE  K  is  available,  it  is  reasonable  to  think  y  / /k  lies  closer 
®  r^+i  c 

to  y  ./vTT.  Hence  a  tangent  line  at  y  . / can  be  used  as  an  approxi- 

•  j  '  *  w  I  ^  '  i  w 

mation  of  g(z)  in  the  neighborhood  of  y  1 / /K  (see  Fig.  2). 


i.e.  ndw  8  *  g’(yr  +1/^T)  (5.1) 

and  a  *  g(yr^+1/»^)  -  (y^^/^)  *  S  , 

and  computing  K,  by  the  same  formula  (2.13)  as  before  will  give  a  new  estimator  K^. 
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The  following  tables  show  the  and  MLE  by  iteration  for  sample  size 


n  *  10  and  n  = 

30,  K  -  3. 

Table  6.  n  * 

10 

*1 

q2 

MLE 

. Error 

.1 

.1 

2.5133 

2.5123 

.0003 

.1 

.2 

2.0379 

2.0369 

.0004 

.1 

.3 

.8669 

.8669 

.0000 

.1 

.4 

.9576 

.9575 

.0000 

.1 

.5 

.9576 

.9575 

.0000 

.2 

.1 

2.S337 

2.5335 

.0000 

.2 

.2 

2. 0605 

2.0603 

.0001 

.2 

.3 

2.0605 

2.0603 

.0001 

.2 

.4 

.8947 

.8944 

.0004 

.2 

.5 

.9915 

.9912 

.0002 

.3 

.1 

2.5335 

2.5335 

.0000 

.3 

.2 

2.0603 

2.0603 

.0000 

.3 

.3 

2.0603 

2.0603 

.0000 

.3 

.4 

.8943 

.8943 

.0000 

.3 

.5 

.9912 

' .9912 

.0000 

.4 

.1 

2.4952 

2.4952 

.0000 

.4 

.2 

2.0163 

2.0163 

.0000 

.4 

.3 

2.0163 

2.0163 

.0000 

V 

.4 

.4 

.8265 

.8264 

.0001 

.4 

.5 

.9117 

.9113 

.0004 

Error  ■  |K£-MLE|/MLE 
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Table  7.  n  ■  30,  K  «  3. 


ql 

<2 

K' 

c 

MLE 

Error 

.1 

.1 

2.94108 

2.94108 

.0 

.1 

.2 

2.39997 

2.39997 

•  .0 

.1 

.3 

2.07875 

2.07874 

.0 

.1 

.4 

1.91790 

1.91788 

.1 

.1 

.5 

1.95178 

1.95175 

.2 

.2 

.1 

2.95687 

2.95687 

.0 

.2 

.2 

2.41740 

2.41738 

.1 

.2 

.3 

2.09744 

2.09742 

.1 

.2 

.4 

1.94050 

1.94047 

.2 

.2 

.5 

1.97899 

1.97894 

.3 

.3 

.1 

2.95151 

2.95148 

.1 

.3 

.2 

2.41096 

2.41096 

.0 

.3 

.3 

2.09010 

2.09010 

.0 

.3 

.4 

1.93130 

1.93127 

.1 

.3 

•5 

1.96801 

1.96796 

.2 

.4 

.1 

2.99032 

2.99031 

.0 

.4 

.2 

2.45209 

2.45207 

.1 

..4 

.3 

2.13266 

2.13261 

.2 

.4 

.4 

1.98169 

1.98160 

.4 

.4 

.5 

2.02922 

2.02906 

.8 

Error  «  (jK^-MLE|/MLE)*(104) 
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Table  8  and  9  show  K*  and  MLE  by  iteration  when  K  *  10.000. 

c 

Table  8.  n  -  10,  K  -  10.0 


ql 

q2 

K' 

c 

MLE 

Error 

.1 

.1 

8.37781 

8.37461 

.00038 

.1 

.2 

6.79320 

6.78998 

.00047 

.1 

.3 

6.79320 

6.78998 

.00047 

.1 

.4 

2.88993 

2.88993 

.00000 

.1 

.5 

3.19218 

3.19195 

.00007 

.2 

.1 

8.44567 

8.44496 

. 00009 

.2 

.2 

6.86851 

6.86776 

. 00011 

.2 

.3 

6.86851 

6.86776 

.00011 

.2 

.4 

2.98254 

2.98119 

.00045 

.2 

.5 

3.30493 

3.30404 

.00027 

.3 

.1 

8.44497 

8.44496 

. 00000 

.3 

.2 

6.86789 

6.86776 

.00002 

.3 

.3 

6.86789 

6.86776 

.00002 

.3 

.4 

2.98122 

2.98119 

.00001 

.3 

.5 

3.30424 

3.30404 

.00006 

.4 

.1 

8.31778 

8.31747 

.00005 

.4 

.2 

6.72130 

6.72129 

.00000 

.4 

.3 

6.72130 

6.72129 

.00000 

.4 

.4 

2.75501 

2.75471 

.00011 

.4 

.5 

3.03910 

3.03787 

.00041 

Error  -  |K^-MLE|/MLE 


Table  9.  n  ■ 

ql 

30,  K. 

q2 

-  10.000 

K* 

c 

MLE 

Error 

.1 

.1 

9.80361 

9.80361 

.00000 

.1 

.2 

7.99990 

7.99989 

. 00000 

.1 

.3 

6.92917 

6.92914 

. 00000 

.1 

.4 

6.39299 

6.39293 

.00001 

.1 

.5 

6.50595 

6.50583 

.00002 

.2 

.1 

9.85625 

9.85623 

. 00000 

.2 

.2 

8.05799 

8.05793 

.00001 

.2 

.3 

6.99148 

6.99141 

. 00001 

.2 

.4 

6.46834 

6.46824 

.00002 

.2 

.5 

6.59664 

6.59646 

.00003 

.3 

.1 

9.83837 

9.83827 

.00001 

.3 

.2 

8.03654 

8.03654 

.00000 

.3 

.3 

6.96701 

6.96699 

.00000 

.3 

.4 

6.43766 

6.43756 

.00001 

.3 

.5 

6.56002 

6.55988 

.00002 

.4 

.1 

9.96773 

9.96770 

.00000 

.4 

.2 

8.17363 

8.17356 

.00001 

.4 

.3 

7.10888 

7.10871 

.00002 

.4 

.4 

6.60563 

6.60535 

.00004 

.4 

.5 

6.76406 

6.76354 

.00008 

Error  -  ( |  K^,-MLE  j /MLE) 
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Conclusion i  The  approximation  1?£  to  an  MLE  which  was  proposed  in  this  sec¬ 
tion  has  some  desirable  advantages  over  the  actual  iteration  procedure  Con¬ 
sidering  the  fact  that  the  iteration  procedure  requires  several  steps,  is 
clearly  time  saving  because  the  procedure  suggested  requires  the  time  equiva¬ 
lent  to  that  which  is  needed  for  a  single  iteration. 

The  iteration  procedure  depends  heavily  on  the  initial  guess  and  as  such 
it  is  the  main  reason  we  do  not  want  to  compare  the  computing  times  required 
by  our  method  and  the  iteration  procedure.  If  the  initial  guess  is  poor  the 
computing  time  required  by  the  iteration  will  be  lengthened.  In  fact,  we 
experienced  several  divergence  of  the  Newton-Raphson' s  method  for  our  simple 
example  though  the  initial  guess  was  chosen  as  2.0  when  the  actual  parameter 
was  3.  is  computed  using  the  previously  obtained  K£  and  hence  it  has  an 
advantage. 

As  regards  the  accuracy  of  the  procedure  suggested,  we  tried  some  more 
simulated  data  obtained  from  the  distribution  with  different  parameters  and 
concluded  that  the  relative  error  to  the  iterated  estimate  is  negligible  in 
most  cases.  (Less  chan  10  for  the  sample  of  size  10  or  more).  Hence  we 
recommend  to  use  when  the  time  saving  is  a  crucial  matter  or  when  there  is 
a  difficulty  of  setting  the  initial  guess. 

At  present  we  are  investigating  the  possibility  of  extending  this  new 
procedure  to  the  results  given  in  the  series  of  papers  by  Tiku. 
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